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Preface 
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Chapter  I:  Introduction 


The  Study  System 

Indian  River  and  Rehobo th  Bay  (Figure  1-1)  are  two  water  bodies  that  form 
part  of  the  interlocked  Delaware  Inland  Bays  system.  Rehobo  th  Bay  is  con¬ 
nected  to  Delaware  Bay  to  the  north  via  a  canal  and  to  Indian  River  Bay  to  die 
south.  Indian  River  Bay  is  connected  to  the  Atlantic  Ocean  on  the  east  via  an 
inlet  and  to  Little  Assawoman  Bay  via  a  canal  to  the  south.  The  western 
portion  of  Indian  River  Bay,  referred  to  as  Indian  River,  terminates  at  the 
Millsboro  dam.  Drainage  area  of  the  system  is  55647  hectares  (Ritter  1986)  of 
which  14339  hectares  is  upstream  of  the  impoundment  at  Millsboro.  The 
basin  contains  one  long-term  stream  gauging  station  (USGS  01484500)  on 
Stockley  Branch.  Mean  flow  for  the  period  of  record  (43  years)  is  0.196  m3 
sec'1  or  1.44  x  KT4  m3  sec'1  hectare'1.  Employing  the  runoff  at  Stockley  to 
characterize  the  remainder  of  the  basin  indicates  a  long-term  basin  mean  flow 
of  8.03  m3  sec'1. 

Surface  area  and  volume  of  the  system  are  7.31  x  107  m2  and  121  x  10s 
m3  respectively.  Mean  depth  is  1.66  m  which  characterizes  most  of  die  system 
except  near  the  inlet  in  which  local  mean  depth  exceeds  ten  meters.  Mean  tide 
range  at  die  inlet  is  1.25  m  (Smullen  1992).  The  tidal  prism  is  51  x  10^  m3 
(Smullen  1992).  The  system  is  well-mixed  from  surface  to  bottom  and  is 
saline  virtually  throughout  its  extern.  Median  salinity  is  22.7  ppt  and  95%  of 
observations  exceed  4.3  ppt.  Lowest  salinities  occur  immediately  downstream 
of  the  Millsboro  dam  during  periods  of  high  runoff.  Residence  time  of  die 
system,  determined  as  volume  divided  by  freshwater  flow  rate  is  lengthy: 

174  days.  An  alternate  estimate  of  residence  time,  volume  divided  by  tidal 
prism  over  the  tidal  period,  is  much  less:  12  days.  Except  near  headwaters 
and  in  constricted  areas  in  which  the  tide  is  dampened,  tidal  flushing  is  more 
effective  than  runoff  in  the  determination  of  volumetric  flows  and  mass  trans¬ 
port  throughout  die  system. 


Objectives 

The  primary  objective  of  this  study  is  to  provide  a  hydrodynamic/water 
quality  model  package  of  the  Indian  River  •  Rehoboth  Bay  system.  The 
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package  is  to  be  suitable  for  development  of  total  maximum  daily  loads,  point- 
source  waste  load  allocations,  and  nonpoint-source  load  allocations  of  nutrients 
and  organic  substances.  Additional  requirements  are  that  the  model  package 
operate  in  a  continuous  multi-year  mode  and  provide  information  on  diurnal 
dissolved  oxygen  variations.  The  model  package  is  also  to  provide  an  orga- 
nized  framework  for  collection  and  employment  of  additional  observations  in 
the  study  system.  The  period  covered  by  the  study  extends  from  January  1988 
to  December  1990. 
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Chapter  II:  Data  Bases 


Hydrographic  Data  Bases 

Bathymetry 

A  hydrographic  survey  of  the  dual  bay  and  inlet  system  was  perforated  in 
1988  by  the  U.S.  Army  Engineer  District,  Philadelphia.  The  survey  was  refer¬ 
enced  to  State  Plane  coordinates  and  all  soundings  were  referenced  to  National 
Geographic  Vertical  Datum  (NGVD).  These  data  were  used  to  determine  the 
depth  in  each  cell  of  the  Indian  River-Rehoboth  Bay  grid  by  applying  a 
3-point  linear  interpolation  scheme  to  the  data  In  addition,  the  data  were  con¬ 
verted  from  State  Plane  coordinates  to  map  inches  for  consistency  between  the 
grid  reference  system  and  bathymetric  reference  system. 


Tide  gages 

Tide  data  referenced  to  NGVD  were  collected  by  the  United  States  Geo¬ 
logic  Survey  at  five  locations  (Figure  2-1).  These  data  (Table  2-1)  were  used 
to  specify  the  inlet  boundary  condition  and  for  calibration  purposes.  Tide 
records  at  die  inlet  were  not  complete  for  the  three-year  model  application 
period,  1988-1990.  Missing  data  were  obtained  by  transformation  of  tide 
records  at  Lewes,  Delaware.  The  transformation  was  accomplished  by  com¬ 
paring  the  inlet  and  Lewes  gages  for  the  time  period  when  data  were  available 
at  both  gages  and  establishing  a  relationship  between  the  two  gages. 


Current  meters 

Velocity  data  employed  in  the  model  calibration  were  collected  at  four 
locations  near  Middle  Island  from  30  June  to  1  July  1988  (Figure  2-1  and 
Table  2-2). 


Figure  2-1 .  Gauge  locations  for  model  calibration 


Wind  Data 

Three  sources  of  wind  data  were  available.  The  State  of  Delaware,  Divi¬ 
sion  of  Soil  and  Water  Conservation,  provided  data  from  an  anemometer 
mounted  above  the  dune  line  at  die  “old”  United  States  Coast  Guard  (USCG) 
station  (Figure  2*1).  The  Delmarva  Power  and  Light  (DPL)  Company  pro¬ 
vided  wind  data  at  the  power  plant  (Figure  2-1).  Weather  data  from  the  Dover 
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Table  2-1 

Tide  Gauge  Data 

Location 

Gauge  Number 

Northing 

Easting 

Time  Periods 

Vines 

(Tl) 

01484549 

202070 

560870 

06-29-88  to  07-05-88 
10-01-88  to  09-30-89 
10-01-89  to  10-13-89 
07-06-90  to  09-30-90 

Pot  nets 

(T2) 

01484605 

222270 

576670 

06-29-88  to  09-30-88 
10-01-88  to  09-07-89 
09-25-89  to  09-30-89 
10-01-89  to  07-17-90 

Massey 

(T3) 

01484680 

227670 

590800 

06-29-88  to  08-15-88 
08-17-88  to  09-30-88 
10-01-88  to  09-30-89 
10-01-89  to  07-17-90 

Dewey 

(T4) 

014846’0 

253070 

594750 

06-29-88  to  07-05-88 
10-01-88  to  09-30-89 
10-31-89  to  09-30-90 
10-01-90  to  11-04-91 

Inlet  (USCG) 

(T5) 

01484683 

222170 

599600 

06-29-88  to  09-30  88 
10-01-88  to  10-09-88 
10-13-88  to  06-19-89 

Table  2-2 

Current  Meter  Data 

Location 

Northing 

Easting 

Time  Periods 

east  of  Middle  Island 
(east  gauge)  (CMl) 

224610 

592800 

06-30-88  to  07-01-88 

east  of  Middle  Island 
(west  gauge)  (CM2) 

224900 

592470 

06-30-88  to  07-01-88 

west  of  Middle  Island 
(east  gauge)  (CM3) 

225600 

591330 

06-30-88  to  07-01-88 

west  of  Middle  Island 
(west  gauge)  (CM4) 

225730 

591200 

06-30-88  to  07-01-88 

Delaware  Air  Force  Base  provided  the  third  source  of  wind  data.  For  the  1988 
calibration,  data  from  the  DPL  anemometer  were  used  because  this  local  data 
source  was  most  complete.  For  the  1989  and  1990  simulations,  weather  data 
from  Dover  was  the  most  complete  data  source  and  was  therefore  utilized. 


Chapter  2  Data  Bases 
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Water  Quality  Data 


STORET 

Data  from  over  60  stations  within  the  Indian  River  -  Rehoboth  Bay  water¬ 
shed  was  retrieved  from  the  STORET  data  base.  Stations  were  plotted  on 
topographic  maps  according  to  reported  longitude  and  latitude.  In  many  cases, 
reported  longitude  and  latitude  were  not  consistent  with  landmarks  in  the  sta¬ 
tion  descriptions.  Longitude  and  latitude  were  then  revised  to  conform  to 
station  descriptions.  Some  stations  could  not  be  accurately  located  and  were 
dropped  as  were  stations  that  contained  no  water  quality  constituents  of  inter¬ 
est.  The  STORET  data  was  divided  into  two  data  sets,  “Upland”  data  from 
freeflowing  streams  and  “Water”  data  collected  in  the  tidal  portions  of  Indian 
River  and  Rehoboth  Bay.  “Upland”  data  was  used  to  compute  distributed 
loads  to  the  system  and  to  characterize  boundary  conditions  in  the  freeflowing 
streams.  “Upland”  data  was  available  for  35  stations  and  was  collected  in  the 
years  1970  -  1991.  All  months  of  the  year  were  represented  in  the  “Upland” 
data  set.  “Water”  observations  were  employed  to  calibrate  and  verify  the 
water  quality  model.  “Water”  data  was  available  for  18  stations  (Figure  2-2) 
and  was  restricted  to  the  years  within  the  study  period,  1988-1990.  Virtually 
all  (97.6%)  of  the  “Water”  data  was  collected  in  the  months  April  -  September. 
The  number  of  observations  at  each  station  ranged  from  2  to  18  and  varied 
according  to  constituent.  STORET  data  employed  in  this  study  is  summarized 
in  Table  2-3. 


University  of  Delaware 

Additional  data  for  calibration  and  verification  of  the  model  was  provided 
by  Dr.  William  Ullman  of  the  College  of  Marine  Studies  (CMS),  University  of 
Delaware.  Protocol  in  the  CMS  study  called  for  sample  collection  at  consis¬ 
tent  salinity  concentrations  rather  than  consistent  physical  locations.  As  a 
consequence  each  observation  was  generally  collected  in  a  unique  location. 
Surveys  were  conducted  in  October  1989  (Figure  2-3)  and  in  March,  May, 
July,  August,  and  September  1990  (Figure  2-4).  The  CMS  data  set  is  also 
summarized  in  Table  2-3. 


Sediment-Water  Fluxes 

Observations  of  sediment-water  fluxes  of  dissolved  oxygen,  ammonium, 
phosphate,  and  nitrate  were  provided  by  Dr.  Sybil  Seitzinger  of  the  Academy 
of  Natural  Sciences.  Observations  were  collected  at  four  stations  (Figure  2-5) 
during  May  and  August  1992. 
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Figure  2-2.  STORET  Data  Locations 


Light  Extinction  Data 

Observations  of  disk  visibility  (secchi  depth)  and  chlorophyll  collected  at 
nine  stations  (Figure  2*6)  during  1985  and  1986  were  obtained  from  an  Aca¬ 
demy  of  Natural  Sciences  Report  (Academy  of  Natural  Sciences,  1988).  These 
observations  were  combined  with  disk  visibility  and  chlorophyll  data  collected 
at  numerous  locations  by  the  College  of  Marine  Studies  from  July  1990 
through  October  1992.  The  combined  data  base  was  employed  in  the  compu¬ 
tation  of  light  extinction. 
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Table  2-3 

Observations  tot  Major  Data  Sets 

Source 

STORET  “Water" 

University  of 

Delaware 

STORET  “Upland" 

Ported 

March  19M  • 
December  1990 

October  1999  - 
September  1990 

1970  - 1991 

Stations 

18 

185 

35 

Temperature 

205 

185 

641 

Salinity 

205 

148 

617 

Ammonium 

205 

185 

623 

Total  Kjetdahl  Nitrogen 

182 

0 

604 

Nitrate 

126 

185 

629 

Phosphate 

78 

185 

271 

Total  Phosphorus 

173 

0 

477 

Dissolved  Oxygen 

205 

0 

631 

BOD5 

51 

0 

481 

Chlorophyll  'a' 

81 

184 

9 

Diurnal  Dissolved  Oxygen  Surveys 

Diurnal  dissolved  oxygen  data  were  obtained  from  two  sources.  The 
Department  of  Natural  Resources  and  Environmental  Control  provided  data 
collected  at  four  locations  (Figure  2-7)  in  July  1991.  Observations  of  tempera¬ 
ture,  salinity  and  dissolved  oxygen  were  conducted  at  fifteen-minute  intervals 
for  24  hours.  Diurnal  data  collected  at  five  stations  (Figure  2-7)  in  August 
1983  were  obtained  from  a  University  of  Delaware  report  (Biggs,  1984). 
Hourly  measures  of  dissolved  oxygen  conducted  for  24  hours  were  supple¬ 
mented  with  less  frequent  measures  of  nitrogen,  phosphorus,  and  chlorophyll. 


Offshore  Water  Quality 

Water  quality  observations  collected  offshore  of  the  Indian  River  inlet  were 
provided  by  the  Region  III  U.S.  Environmental  Protection  Agency.  Data  were 
from  May  -  November,  1988  •  1990,  and  were  collected  as  part  of  the  EPA 
Coastal  Eutrophication  Surveys.  The  observations  were  employed  in  the  speci¬ 
fication  of  oceanic  boundary  conditions  for  temperature,  salinity,  ammonium, 
nitrate,  phosphate,  chlorophyll,  and  dissolved  oxygen. 
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Figure  2-3.  University  of  Delaware  Data  Stations,  1989 
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Figure  2-4.  University  of  Delaware  Data  Stations,  1990 
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Figure  2-6.  Location  of  Academy  of  Natural  Sciences  Disk  Visbility  Measures 


Figure  2-7.  Location  of  Diurnal  Dissolved  Oxygen  Surveys 
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Chapter  III:  Flows  and  Loads 


Flows 

Freshwater  Runoff 

The  volumetric  freshwater  runoff  rate  was  one  of  the  forcing  functions 
required  to  drive  the  hydrodynamic  portion  of  the  Indian  River/Rehoboth  Bay 
model.  Both  the  time  series  and  spatial  distribution  of  runoff  were  necessary. 
Runoff  volume  was  also  required  to  compute  the  distributed  loads  of  nutrients 
and  organic  matter  to  die  system.  The  drainage  basin  contained  only  one 
stream  gauge,  on  Stockley  Branch  (Figure  3-1),  active  during  the  entire  study 
period.  The  strategy  adopted  to  get  flow  throughout  the  basin  was  to  convert 
Stockley  volumetric  flow  (e.g.  nr  sec'1)  to  flow  per  unit  area  (e.g.  cm  day'1). 
This  flow  per  unit  area  was  multiplied  by  drainage  areas  of  subbasins  within 
the  watershed  to  get  the  volumetric  flow  in  each  subbasin.  Subbasin  locations 
and  areas  were  obtained  from  Ritter  (1986)  who  identified  16  subbasins  in  the 
Indian  River/Rehoboth  Bay  watershed.  Of  these,  three  contributed  to  Mills- 
boro  Pond  upstream  of  die  dam  and  were  combined  by  us  into  one.  Land 
areas  draining  directly  into  Indian  River  and  Rehoboth  Bay  were  named  as 
individual  watersheds  by  Ritter.  We  allocated  the  direct  discharge  area  to 
tributary  subbasins  so  that  the  total  number  of  subbasins  considered  in  this 
study  was  twelve  (Table  3-1).  Flows  from  each  subbasin  were  input  to  the 
model  at  the  discrete  location  (Figure  3-1)  at  which  the  subbasin  tributary 
entered  the  receiving  water.  Flows  were  updated  on  a  daily  basis  throughout 
the  three-year  study  period. 


Millsboro  Pond 

Flow  from  the  largest  of  the  subbasins  (Figure  3-2)  altered  the  tidal  system 
across  the  Millsboro  spillway.  A  gauging  station  (USGS  0148453525) 
operated  at  die  spillway  during  a  portion  of  the  study.  Comparison  of  flows  at 
Stockley  and  Millsboro  (on  a  per  unit  area  basis)  showed  die  impoundment 
exerted  a  smoothing  effect  (Figure  3-3).  Peak  flows  at  Millsboro  were  damp¬ 
ened  relative  to  die  Stockley  gauge  while  minimum  flow  from  die  impound¬ 
ment  was  higher  than  gauged  at  the  upland  station.  Due  to  the  smoothing 
effect,  the  Stockley  flow  per  unit  area  was  not  used  to  provide  flows  at 
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Figure  3-1 .  Locations  of  distributed  flows 


Millsboro.  Rather,  a  regression  relationship  was  derived  that  related  flow  at 
Millsboio  to  flow  cm  the  Nanticoke  River  near  Bridgeville.  The  Nanticoke 
gauge  (USGS  01487000)  was  nearby  but  outside  die  Indian  River  watershed. 
The  Nanticoke  drainage  area  (19529  hectares)  was  an  order  of  magnitude 
larger  than  the  Stockley  drainage  area  (1357  hectares).  The  large  drainage  area 
of  the  Nanticoke  gauge  apparently  exerted  a  smoothing  effect  on  flow  similar 
to  the  smoothing  effect  of  the  Millsboro  impoundment.  The  relationship  was: 

FMB  *  1.22  FNT  -  0.133  t3"1) 

in  which: 

FMB  =  flow  at  Millsboro  (cm) 

FNT  =  flow  at  Nanticoke  (cm) 
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Figure  3-3.  Comparison  of  flow  per  unit  area  at  Stockley  and  MiHsboro 


The  R2  for  the  relationship  was  0.95  for  monthly-average  flows. 


Delmarva  Power  and  Light  Diversion 

The  Delmarva  Power  and  Light  Company  operates  a  power  plant  on  the 
southern  shore  of  Indian  River.  Three  once-through  cooling  units  divert  water 
from  Indian  River,  through  the  power  plant  condensers,  and  into  Island  Creek. 
Discharge  water  flows  down  die  creek  and  rejoins  Indian  River  downstream  of 
the  intake,  close  to  the  location  where  the  riverine  portion  opens  out  into 
Indian  River  Bay  (Figure  3-1).  The  cooling  water  diversion  was  included  in 
the  hydrodynamic  model.  Flow  through  the  power  plant,  at  monthly  intervals, 
was  obtained  from  discharge  monitoring  records  (DMRs)  provided  by  the 
Department  of  Natural  Resources  and  Economic  Conservation.  The  flow 
diversion  was  12  to  16  m3  sec1  (Figure  3-4).  Interpolated  flows  were 
employed  to  fill  gaps  in  the  DMRs. 


Relation  to  Groundwater  Row 

Roughly  80%  of  stream  flow  in  the  Indian  River/Rehoboth  Bay  watershed 
is  base  flow  contribution  from  groundwater  (Johnston  1976).  The  runoff  per 
unit  area  records  derived  from  the  Stockley  gauge  included  both  groundwater 
and  overland  runoff.  Extension  of  the  Stockley  record  to  the  entire  watershed 
extended  both  overland  and  groundwater  runoff  volumes.  Independent  esti¬ 
mates  (Andres  1992)  have  been  made  of  the  groundwater  contribution  to 
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Figure  3*4.  Volumetric  flow  diverted  through  Delmarva  Power  and  Light  cooling  units 


Indian  River/Rehoboth  Bay  from  drainage  immediately  adjacent  to  the  shore¬ 
line.  Since  this  area  was  represented  in  our  twelve  subbasins,  no  additional 
accounting  for  groundwater  was  required. 


Flows  Through  Navigation  Canals 

Canals  connect  Rehobo th  Bay  with  Delaware  Bay  and  Indian  River  Bay 
with  Little  Assawoman  Bay.  At  present,  no  information  exists  on  flow  and 
material  exchange  through  these  canals.  As  a  first  approximation,  we  assumed 
no  net  flow  occurs  through  the  canals. 


Hydrologic  Characterization 

Runoff  at  the  Stockley  gauge  during  the  study  period  is  summarized  in 
Figure  3-5.  Two  years,  1988  and  1990,  had  typical  hydrographs.  Flow  was 
highest  during  the  first  half  of  the  year  and  lowest  during  spring  and  fall.  The 
hydrograph  for  1989  was  unusual,  however.  Peak  flow  occurred  during 
August  and  flow  in  all  months  exceeded  the  long-term  average.  Recurrence 
relationships  (Table  3-2)  indicate  the  flows  that  occurred  in  1989  are  exceeded 
by  less  than  8%  of  the  years  on  record.  Flows  that  occurred  during  summer 
1989  are  exceeded  by  less  than  2%  of  the  years  on  record.  Calendar  year 
1989  was  extremely  wet  By  contrast,  annual  flow  in  1990  was  close  to  long¬ 
term  median  flow  while  1988  annual  flow  was  below  average.  Summer  flows 


Chapter  3  Flows  and  Loads 


3-5 


45 

4 

35 

3 


>-  25 

s 

8  2 


15 


1 


05 

0 


SSSSS8SR2SKS8SSSKS82K 

OO  CO  CO  CO  CO  oo  CO  CO  OO  OO  CO  CO 
I  I  I  I  I  I  I  t  I  I  I  I 


<n  txt  cn  cr>  O)  cr>  01  0>  cr>  01  01  01  o  o  o  o  o  o  o  o  o  o  o  o 
aOCOaOOOQOOOaOQOflOCOOOOO0>0)0)0)O)0)0)0)0)0)0l0> 
I  I  I  I  I  I  I  I  I  I 


-15  s 


1  o'  a-r 
ra  a» 


Figure  3-5.  Monthly  flow  per  unit  area  at  Stockley  1988  - 1990 


Table  3-2 

Hydrologic  Characterization 

rtrlOG 

Percent  of  Years  That  Exceed  This  One* 

Characterization 

Twelve  Months,  1988 

70  2 

Dry 

Twelve  Months,  1989 

7.4 

Wet 

Twelve  Months,  1990 

54  2 

Average 

July  -  September,  1988 

43.6 

Average 

July  -  September,  1989 

1.1 

Wet 

July  -  September,  1990 

39.3 

Average 

||  1  Obtained  from  SMullen  1992 

in  1988  and  1990  were  similar  and  close  enough  to  the  median  to  be  character¬ 
ized  as  "average”. 
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Loads 


The  water  quality  model  requires  nutrient  and  organic  loads  to  Indian 
River/Rehoboth  Bay  as  forcing  functions  for  computation  of  receiving  water 
quality.  Loads  can  be  divided  into  three  classes:  Distributed  Loads,  Point- 
Source  Loads,  and  Atmospheric  Loads.  Distributed  loads  enter  the  system  as 
runoff  from  the  watersheds  that  drain  into  the  tidal  waters.  Point-source  loads 
originate  in  industries  and  treatment  plants  along  the  shoreline.  Atmospheric 
loads  are  deposited  directly  on  the  water  surface  in  rainfall  and  as  dryfall. 
Both  the  time  series  and  locations  of  these  loads  are  required  for  the  model. 


Distributed  Loads 

Distributed  loads  were  computed  as  the  product  of  runoff  volume  and  con¬ 
centration  of  nutrient  or  organic  substance.  Runoff  volumes  were  the  same  as 
those  used  to  drive  tire  hydrodynamic  model.  Monthly  loads  were  computed 
for  each  subbasin  and  input  to  the  model  at  the  same  locations  as  the  flows 
(Figure  3-1).  Concentrations  in  the  runoff  were  obtained  from  analysis  of  the 
“Upland”  stations  in  the  STORET  data  base.  Separate  analyses  were  con¬ 
ducted  for  the  Millsboro  subbasin  and  for  the  remaining  watersheds.  We 
conducted  separate  analyses  to  allow  for  concentration  differences  between 
freeflowing  streams  and  water  leaving  the  impoundment. 

Freeflowing  Streams.  Our  analysis  of  the  data  indicated  no  relationship 
between  concentration  and  location  or  flow.  Observations  of  most  substances 
exhibited  a  high  degree  of  variability.  Mean  concentrations  were  unduly  influ¬ 
enced  by  a  few  extreme  observations.  We  characterized  concentrations  in  the 
runoff  (Table  3-3)  as  the  median  value  of  all  observations  since  the  median 
statistic  is  less  influenced  by  extreme  values  than  the  mean.  No  temporal  or 
spatial  variability  in  concentration  was  considered  except  for  nitrate  which 
exhibited  a  temporal  pattern;  concentration  was  generally  higher  in  late  autumn 
and  winter  than  in  summer  (Figure  3-6).  Monthly  median  concentrations  were 
employed  to  compute  nitrate  loads. 

Millsboro  Spillway.  Analysis  of  concentrations  in  water  leaving  the  Mills¬ 
boro  impoundment  parallelled  analysis  of  concentration  in  freeflowing  streams. 
No  relationship  of  concentration  to  flow  or  season  was  evident  except  for 
nitrate  which  showed  a  temporal  pattern  similar  to  freeflowing  streams  (Fig¬ 
ure  3-6).  Concentrations  were  characterized  as  the  median  of  all  observations 
(Table  3-4)  except  for  nitrate.  For  nitrate,  monthly  median  values  were 
employed. 

Total  Organic  Carbon.  The  water  quality  model  employs  organic  carbon 
as  a  state  variable.  Runoff  observations  included  BODS,  but  not  organic  car¬ 
bon.  No  direct  conversion  of  BOD  to  organic  carbon  exists  but  empirical 
relationships  can  sometimes  be  found.  We  examined  a  data  base  of  total 
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Table  3-3 

Concentrations  In  Freeflowing  Streams 

Ammonium 

Nitrate 

Total 

Kjoldahl 

Nitrogen 

Phosphate 

Organic 

Phosphorus 

Total 

Organic 

Carbon 

Month 

gmm4 

pm4 

gmm4 

gmm4 

gmm4 

gmm4 

Jan 

0.1 

2.98 

0.88 

0.11 

0.03 

9 

Fab 

0.1 

2.58 

0.88 

0.11 

0.03 

9 

Mar 

0.1 

2.65 

0.88 

o.n 

0.03 

9 

Apr 

0.1 

2.04 

0.88 

o.n 

0.03 

9 

May 

0.1 

1.62 

0.88 

0.11 

0.03 

9 

Jun 

0.1 

0.86 

0.88 

O.n 

0.03 

9 

Jul 

0.1 

0.23 

0.88 

0.11 

0.03 

9 

Aug 

0.1 

1.13 

0.88 

0.11 

0.03 

9 

Sap 

0.1 

1.69 

0.88 

0.11 

0.03 

9 

Oct 

0.1 

0.39 

0.88 

0.11 

0.03 

9 

Nov 

0.1 

1 

0.88 

0.11 

0.03 

9 

Dec 

0.1 

2.43 

0.88 

0.11 

0.03 

9 

organic  carton  (TOC)  and  CBOD5  observations  collected  in  freeflowing 
streams  tributary  to  the  upper  Potomac  River.  A  loose  curvilinear  relationship 
was  evident  (Figure  3-7).  TOC  concentration  at  Millsboro  and  in  freeflowing 
streams  was  obtained  visually  from  the  figure  as  the  TOC  concentration 
(9  gm  m'3)  corresponding  to  median  BOD5  concentration  (2.5  gm  m'3)  in  die 
“Upland"  data  set. 


Point-source  Loads 

Thirteen  point  sources  exist  within  the  Indian  River/Rehoboth  Beach  water¬ 
shed  (Figure  3-8).  One  of  these,  the  Lewes  STP,  discharges  near  the  northern 
terminus  of  the  canal  that  connects  Delaware  and  Rehoboth  Bays.  Consistent 
with  our  assumption  of  no  net  flow  through  the  canal,  we  assumed  that  Lewes 
effluent  does  not  travel  the  length  of  the  canal  to  Rehoboth  Bay.  The 
Rehoboth  Beach  WWTF  and  lesser  point  sources  that  discharge  near  the  south¬ 
ern  terminus  of  the  canal  were  included  as  Rehoboth  Bay  point  sources.  The 
Georgetown  STP  discharges  into  a  tributary  of  the  Millsboro  impoundment. 
Loads  from  Georgetown  do  not  require  explicit  treatment  since  they  are 
included  in  the  distributed  load  assigned  to  the  Millsboro  spillway.  (Analysis 
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Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec 


Figure  3-6.  Monthly  median  nitrate  concentration  at  Millsboro  spillway  and  in  freeflowing 
streams 


indicates  Georgetown  loads  are  less  than  5%  of  the  total  load  at  Millsboro.) 
Deletion  of  Lewes  and  Georgetown  leaves  11  point  sources  for  consideration. 

Flows  and  Concentrations.  Loads  from  each  point  source  were  computed 
as  the  product  of  volumetric  flow  rate  and  substance  concentration.  Flow  and 
concentration  data  were  supplied  by  DNREC.  Flow  data  was  mostly  complete, 
on  a  monthly  basis,  especially  for  die  larger  point  sources.  Flows  in  missing 
months  were  obtained  by  substituting  flows  from  die  same  months  in  alternate 
years.  Concentration  observations  were  sporadic,  however.  For  some  sources, 
data  were  available  only  for  1991  or  1992,  after  the  study  period.  We  found 
no  rational  basis  for  assigning  temporal  variability  to  the  point-source  concen¬ 
trations.  Mean  concentration  of  all  available  observations  was  assigned  to  each 
point  source  (Table  3-5).  Consequently,  temporal  variations  in  computed 
point-source  loads  were  due  exclusively  to  variations  in  flow  rather  than 
concentration. 

Total  Organic  Carbon.  The  water  quality  model  employs  organic  carbon 
as  a  state  variable.  Point-source  monitoring  included  BOD5,  but  not  organic 
carbon.  Following  the  pattern  set  for  distributed  loads,  we  attempted  to  find 
an  empirical  relationship  between  die  two  substances.  We  examined  a  data 
base  of  total  organic  carbon  (TOC)  and  CBOD5  observations  collected  in  point 
sources  discharging  to  the  upper  Potomac  River.  No  relationship  was  evident 
(Figure  3-9).  We  selected  from  the  data  base  the  median  concentration,  18  gm 
m'3,  to  characterize  effluent  TOC  in  Indian  River/Rehoboth  Bay. 
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Table  3*4 

Concentrations  In  Millaboro  Spillway 


Ammonium 

Nitrate 

Total  KjsMahl 
Nitrogen 

Phoaphata 

-i 

Organic 

Phosphorus 

Total 

Organic 

Carbon 

gmm4 

pm4 

gm  m4 

gmm4 

gmm4 

gmm4 

0.1 

3.07 

0.8 

0.03 

0.02 

9 

0.1 

3.03 

0.8 

0.03 

0.02 

9 

0.1 

2.66 

0.8 

0.03 

0.02 

9 

0.1 

2.23 

0.8 

0.03 

0.02 

9 

0.1 

2.45 

0.8 

0.03 

0.02 

9 

0.1 

2.13 

0.8 

0.03 

0.02 

9 

0.1 

1.35 

0.8 

0.03 

0.02 

9 

0.1 

0.98 

0.8 

0.03 

0.02 

9 

0.1 

1.7 

0.8 

0.03 

0.02 

9 

0.1 

1.76 

0.8 

0.03 

0.02 

9 

0.1 

1.76 

0.8 

0.03 

0.02 

9 

0.1 

25 

0.8 

0.03 

0.02 

9 

Delmarva  Power  and  Light  Information  and  data  on  Delmarva  Power 
and  Light  operations  were  supplied  by  Mr.  Robert  Jubic  of  DP&L.  DP&L  has 
several  operations  which  require  consideration  in  computation  of  point-source 
loads.  The  largest  flow  through  the  plant  is  cooling  water  for  Units  1-3.  No 
material  is  added  to  the  once-through  cooling  water.  We  assumed,  however, 
that  viable  algae  do  not  survive  die  trip  through  the  condensers.  Code  was 
installed  in  the  model  so  that  algae  withdrawn  through  the  diversion  from 
Indian  River  were  discharged  as  equivalent  amounts  of  nutrients  and  organic 
substances  into  Island  Creek.  No  material  was  added  to  the  cooling  water, 
however. 

Water  increases  in  temperature  as  it  passes  through  the  condensers  so  that 
discharge  temperature  is  higher  than  withdrawal  temperature.  Cooling  of  dis¬ 
charge  water  occurs  in  Island  Creek  before  the  water  rejoins  Indian  River. 
Since  our  model  does  not  include  Island  Creek  we  could  not  use  plant  records 
of  discharge  temperature.  We  had  no  way  to  compute  the  cooling  effect  in  the 
Creek.  Instead  we  used  STORET  data,  supplemented  by  University  of  Dela¬ 
ware  observations,  to  derive  temperature  differences  between  observations 
collected  near  the  DP&L  intake  and  die  mouth  of  Island  Creek.  Temperature 
increase  between  the  intake  and  creek  mouth  ranged  from  0.95  to  7.5  C°  (Fig¬ 
ure  3-10).  Average  increase  was  4.2  C°.  Code  was  installed  in  the  model  to 
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TOC  VS  BOD  IN  FREEFLOWING  STREAMS 


Figure  3-7.  Total  organic  carbon  versus  CBOD5  in  freef lowing  streams  tributary 
to  the  Potomac  River  below  Washington,  DC 


increase  by  the  mean  value  temperature  of  water  diverted  through  the  DP&L 
plant 

A  small  portion  of  condenser  discharge  is  diverted  into  the  Unit  4  cooling 
tower.  The  tower  operates  through  evaporative  cooling  so  that  substance  con¬ 
centrations  increase  in  the  cooling  water.  When  solids  concentration  in  die 
water  increases  by  a  factor  of  three  to  four,  die  tower  is  "blown  down” 
(emptied)  via  Discharge  027.  Total  nitrogen  concentration  in  “Mow  down” 
water  is  three  to  four  times  higher  than  intake  water  (Figure  3-1 1),  consistent 
with  the  increase  in  solids.  Nitrate  concentration  increases  by  an  order  of 
magnitude,  however,  and  total  Kjeldahl  nitrogen  concentration  decreases 
despite  die  evaporative  water  loss.  The  data  suggest  mineralization  of  organic 
nitrogen  and  subsequent  nitrification  occur  in  the  tower.  To  account  for  this 
effect  we  converted  appropriate  amounts  of  organic  nitrogen  in  DP&L  diver¬ 
sion  water  (-30  kg  day'1)  to  nitrate.  No  net  addition  or  subtraction  of  total 
nitrogen  was  considered,  however. 

For  two  of  die  three  study  years,  1988  and  1989,  phosphorus  detergent  was 
added  to  the  cooling  tower  water.  Moodily  data  on  phosphorus  use  was  pro¬ 
vided  by  DP&L.  We  added  these  phosphorus  loads  to  the  DP&L  diversion 
water.  Phosphorus  use  at  the  plant  ceased  after  December,  1989. 

No  data  was  available  to  assess  the  effect  of  die  cooling  tower  mi  total 
organic  carbon.  We  assumed  no  net  addition  or  change  through  the  tower. 
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Figure  3-8.  Point  sources  in  the  Indian  River/Rehoboth  Bay  drainage  basin 


The  DP&L  plant  has  a  small  sanitary  facility  (Discharge  028).  This  facility 
was  considered  separately  from  die  cooling  water.  Analysis  and  treatment  of 
the  sanitary  facility  was  identical  to  the  other  point  sources  in  the  system. 


Atmospheric  Loads 

Data  to  compute  atmospheric  nutrient  loads,  collected  at  nearby  Cape 
Henlopen  Delaware,  were  supplied  by  Dr.  Joseph  Scudlark  of  die  University  of 
Delaware.  Mean  ammonium  load  in  rainfall  was  2.6  kg  N  hectare'1  year1. 
Mean  nitrate  load  was  14.1  kg  N  hectare'1  year1.  Organic  nitrogen  load  was 
estimated  as  15%  of  the  inorganic  load.  No  atmospheric  load  of  phosphorus 
was  detected.  A  “rule  of  thumb"  in  consideration  of  atmospheric  loads  is  that 
dryfall,  difficult  or  impossible  to  measure,  equals  wetfalL  To  account  for 
dryfall,  the  measures  and  estimates  of  nitrogen  in  rainfall  were  doubled. 
Atmospheric  loads  were  applied  uniformly  throughout  the  study  period  to  all 
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Table  3*5  J 

Point-Source  Concentrations  j 

Ammonium 

se. —  a- 

Nttrtt* 

Total  KJsktahl 
Nitrogen 

Phosphate 

Organic 

rTIwSp. 

Total 

Org. 

Carh. 

gmm4 

gmm4 

gmm4 

gmm4 

gmm4 

gmm4 

Delaware 
Seashore 
State  Park 

7  3 

17.3 

12 

3.8 

0.6 

18 

Frankford 

Elementary 

School 

9.86 

22.7 

9.86 

53 

1 

18 

Town  of 
Millsboro 

7.1 

2.6 

8.7 

2.4 

0.3 

18 

Townsend’s 

Inc. 

1 

21.8 

6 

0.1 

0 

18 

Vlassic 

Food 

0.06 

3.1 

1.36 

0.064 

0.066 

18 

Colonial 

East  Mobile 
Home  Pk. 

9.7 

6.85 

9.7 

22 

0.4 

18 

Rehoboth 

Beach 

WWTF 

0.95 

0.62 

5.76 

3.46 

0.61 

18 

Delaware 

State 

Housing 

Authority 

1.15 

28.8 

1.15 

435 

0.75 

18 

Bayshore 

Mobile 

Home  Pk. 

73 

173 

12 

3.8 

0.6 

18 

Colonial 

Estates 

7.9 

10.7 

7.9 

3.1 

03 

18 

Deimarva 
Power  and 
Light  (028) 

0.9 

153 

3.1 

136 

0.14 

18 

portions  of  the  water  surface.  No  attempt  was  made  to  account  for  spatial  or 
temporal  variations  in  load.  The  total  atmospheric  nitrogen  load  to  the  water 
surface  was  765  kg  day'1. 
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TOC  VS  BOD  IN  POINT  SOURCES 


Figure  3-9.  Total  organic  carbon  versus  CBOD5  in  point  sources  discharging 
to  the  Potomac  River  below  Washington  DC 


Summary  of  Loads 

Distributed  Loads.  A  summary  of  the  subbasin  loads,  averaged  through¬ 
out  the  study  period  (Table  3-6,  Figure  3-12),  indicates  that  nitrogen  loads  are 
distributed  roughly  in  accordance  with  subbasin  area.  Millsboro  is  die  largest 
subbasin  and  die  largest  distributed  nitrogen  source.  Millsboro  is  also  the 
largest  distributed  phosphorus  source  (Figure  3-13)  but  is  not  as  dominant  as 
for  nitrogen.  The  relatively  large  distributed  phosphorus  loads  from  die  free- 
flowing  subbasins  occur  because  total  phosphorus  concentration  attributed  to 
these  basins  is  nearly  three  times  huger  than  concentration  at  Millsboro.  The 
concentration  disparity  suggests  phosphorus  settling  occurs  in  the  pond  but  the 
apparent  differences  may  also  be  an  artifact  of  the  data  analysis.  Alternate 
treatments  of  the  data  can  be  conducted  that  indicate  no  difference  in  phospho¬ 
rus  concentration.  Total  organic  carbon  loads  from  die  subbasins  (Figure  3-14) 
reflect  exactly  die  relative  areas  of  the  basins  since  concentration  is  uniform 
across  all  subbasins.  Millsboro  is  the  largest  distributed  source  of  total  organic 
carbon. 

Point-Source  Loads.  A  summary  of  the  point-source  loads,  averaged 
throughout  the  study  period  (Table  3-7),  indicates  that  the  Townsend’s  plant  is 
by  far  the  largest  point  source  of  nitrogen  (Figure  3-13).  Rehoboth  Beach 
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Figure  3-1 1 .  Nitrogen  concentrations  in  cooling  tower  intake  and  discharge 
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Table  3-6 

Main  Distributed  Loads  1988-1990 

Total  NHragsn 

ToW  Phosphorus 

TottlOrganie  § 
Carton  1 

Basin 

kg  dsf' 

kg  day* 

kfltoy  I 

Lawes-Rehoboth  Canal 

96 

5.16 

33, 

Lova  Creak 

155 

8.34 

536 

Hatting  Creak 

177 

9.53 

613 

Guinea  Creak 

88 

4.76 

306 

Lingo  Creak 

39 

2.13 

137 

Swan  Creek 

155 

837 

538 

Iron  Branch 

170 

9.15 

588 

Pepper  Creak 

113 

6.07 

390 

Vines  Creak 

109 

5.85 

376 

Blackwater  Creak 

94 

3.05 

325 

White  Creek 

89 

4.78 

307 

MiHsboro  Spillway 

717 

11.90 

2143 

Total 

2000 

81.09 

6591  | 

WWTF  dominates  the  phosphorus  point  sources  (Figure  3-16)  and  is  the 
second  largest  nitrogen  contributor.  Town  of  Millsboro  is  the  second  largest 
phosphorus  contributor  followed  by  the  DP&L  “blow  down"  water  (027). 

Note  that  the  loads  for  DP&L  are  averaged  across  the  three-year  study  period. 
For  two  years,  1988  and  1989,  the  DP&L  load  was  roughly  equivalent  to  die 
Town  of  Millsboro.  The  DP&L  toad  in  1990  was  zero,  however,  and  does  not 
currently  exist  Townsend’s  and  Rehoboth  Beach  also  lead  in  point-source 
total  organic  carbon  (Figure  3-17). 

Relative  Loading.  Distributed  loads  comprise  die  largest  nitrogen  source 
to  the  system  (Figure  3-18).  In  1988  and  1990,  largest  nitrogen  loads  occurred 
in  the  first  six  months.  This  pattern  reflected  both  the  hydrograph  and  the 
temporal  pattern  in  nitrate  concentration.  In  1989,  the  extremely  wet  summer 
was  mirrored  by  unusually  large  nitrogen  loads  in  the  summer  months.  Atmo¬ 
spheric  nitrogen  loads  are  the  next  largest  source.  During  the  dry  summer 
months  of  1988  and  1990,  periods  of  high  algal  and  low  dissolved  oxygen 
concentration,  atmospheric  toads  equaled  or  exceeded  distributed  loads.  Point 
sources  are  die  least  of  the  nitrogen  loads  although  point-source  toads 
approached  distributed  loads  in  magnitude  during  the  dry  summer  of  1988. 
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Figure  3-13.  Mean  total  phosphorus  load  from  subbasins 
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2500.00 


Figure  3-14.  Mean  total  organic  carbon  load  from  subbasins 


Distributed  loads  are  also  the  largest  source  of  phosphorus  (Figure  3-19). 
The  distributed  phosphorus  load  minors  exactly  the  hydrograph  since  no  tem¬ 
poral  variation  in  phosphorus  concentration  is  considered.  In  1988  and  1990, 
largest  phosphorus  loads  occurred  in  the  first  six  months  of  the  year.  In  1989, 
largest  phosphorus  loads  occurred  in  summer  due  to  die  unusually  large  sum¬ 
mer  flows.  Point-source  phosphorus  loads  are  less  than  distributed  loads 
although  during  the  dry  summer  of  1988  point-source  loads  exceeded  distri¬ 
buted  loads.  The  point-source  loads  exhibit  a  periodicity.  Highest  loads  occur 
in  June,  July,  and  August  due  to  flow  variations  at  the  Rehoboth  Beach  facil¬ 
ity,  die  largest  phosphorus  point  source.  Flows  in  summer  months  are  three  or 
four  times  greater  than  in  winter  months. 


As  with  nitrogen  and  phosphorus,  distributed  total  organic  carbon  loads 
exceed  point-source  loads  (Figure  3-20).  The  distributed  loads  dominate 
throughout  the  year  in  dry  or  wet  hydrology. 

Comparison  to  Previous  Estimates.  "Desktop”  estimates  of  loads  to 
Indian  River/Rehoboth  Bay  for  “Wet”,  "Dry”,  and  “Normal”  years  were  com¬ 
pleted  by  Ritter  (1986).  We  compared  our  loads  to  the  previous  estimates. 

Our  loads  for  1988,  1989,  and  1990  were  compared  to  previous  estimates  for 
“Dry”,  “Wet”,  and  “Normal"  years  respectively  (Table  3-8).  Our  total  nitrogen 
loads  are  comparable  to  previous  estimates  but  average  roughly  13%  less. 
Comparison  of  loads  by  category  indicates  rough  equivalence  of  distributed 
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Table  3-7 

Mean  Point-Source  Loads  1988-1990 

Total  Mtrogan 

Total 

Phosphorus 

Total  Organic 
Cartoon 

Source 

kg  day 

kg  day 

kg  day 

Delaware  Seashore  State  Park 

1.37 

0.30 

1.23 

Delmarva  Power  and  Light  (027) 

0.00 

1.64 

0.00 

Frankford  Elementary  School 

0.27 

0.05 

0.15 

Town  of  Millsboro 

11.87 

2.84 

18.91 

Townsend’s  Inc. 

145.23 

0.52 

94.04 

Vlassic  Food 

4.79 

0.14 

19.34 

Colonial  East  Mobile  Home  Pk. 

1.14 

0.18 

1.25 

Rehoboth  Beach  WWTF 

25.76 

16.43 

72.67 

Delaware  State  Housing 
Authority 

0.94 

0.16 

057 

Bayshore  Mobile  Home  Pk. 

0.26 

0.06 

0.24 

Colonial  Estates 

0.90 

0.17 

0.87 

Delmarva  Power  and  Light  (028) 

0.53 

0.04 

052 

Georgetown' 

34.37 

0.48 

27.18 

Lewes  STP* 

32.29 

7.77 

3750 

Total  Contributing  Point  Sources 

193.06 

22.54 

209.77 

1  Load  from  this  point  source  is  included  in  the  distributed  bad  from 
subbasin. 

2  Not  considered  as  a  load  to  Rehoboth  Bay.  Shown  only  for  compa 

the  Millsboro 

rison. 

nitrogen  loads.  Comparison  of  other  categories  indicates  significant  differ¬ 
ences,  however.  Our  point-source  nitrogen  loads  are  only  a  third  of  previous 
estimates.  This  difference  may  be  due  to  a  reduction  in  point-source  loads 
since  the  previous  estimates  were  completed  or  else  improved  quantification 
due  to  availability  of  additional,  contemporary  data.  Our  atmospheric  nitrogen 
loads  are  triple  the  previous  estimates.  In  view  of  the  nature  and  extent  of  the 
Cape  Henlopen  data  base,  current  estimates  must  be  regarded  as  the  best  avail¬ 
able.  The  differences  in  current  and  previous  estimates  of  point-source  and 
atmospheric  loads  offset  Total  nitrogen  loads  would  be  nearly  equivalent 
except  for  additional  loads  considered  by  Ritter.  The  excess  in  previous  esti¬ 
mates  over  current  estimates  is  due  to  septic  tank  loads  considered  by  Ritter 
but  omitted  from  this  study. 
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Figure  3-1 6.  Mean  total  phosphorus  load  from  point  sources 
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Figure  3-1 7.  Mean  t<;r  ,  nanic  carton  toad  from  point  sources 


Figure  3-18.  Time  series  of  point-source,  distributed,  and  atmospheric  nitrogen  loads 
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Figure  3*19.  Time  series  of  point-source  and  distributed  phosphorus  loads 


Figure  3-20.  Time  series  of  point-source  and  distributed  total  organic  carbon  loads 
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Table  3-8 

Comparison  of  Loads 

Nitrogen 

IMS 
kg  day’ 

Dry 

kg  day' 

1M9 
kg  day' 

— 

Wet 

kg  day’ 

1M0 
kg  day' 

Normal 
kg  day’ 

Point 

Sources 

209 

634 

204 

634 

166 

634 

Distributed 

Loads 

1300 

1091 

2683 

3240 

2017 

2186 

Atmospheric 

765 

171 

765 

342 

765 

174 

Other 

0 

515 

0 

515 

0 

515 

Total 

2275 

2410 

3652 

4730 

2948 

3509 

Phosphorus 

Dry 

Wet 

Normal 

Point 

Sources 

24 

72 

23 

72 

21 

72 

Distributed 

Loads 

48 

52 

121 

159 

74 

104 

Atmospheric 

0 

11 

0 

21 

0 

16 

Other 

0 

14 

0 

14 

0 

14 

Total 

72 

149 

144 

267 

96 

206 

Current  estimates  of  total  phosphorus  loads  are  less  than  previous  estimates 
for  all  loading  categories.  The  sum  of  all  current  loads  is  half  the  previous 
estimate.  Differences  in  distributed  loads  are  small  relative  to  die  uncertainty 
in  the  load  estimates.  Differences  in  point-source  and  atmospheric  loading 
likely  occur  due  to  the  availability  to  us  of  additional,  contemporary  data.  As 
with  nitrogen,  phosphorus  loads  from  sources  considered  by  Ritter,  primarily 
septic  tanks,  were  omitted  in  this  study. 
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Chapter  IV:  Description  of  the 
Hydrodynamic  Model 


The  numerical  hydrodynamic  model  CH3D  (Curvilinear  Hydrodynamics  in 
Three  Dimensions)  was  used  to  provide  detailed  hydrodynamic  flow  field 
information  as  input  to  die  water  quality  model.  The  basic  model  was  devel¬ 
oped  by  Sheng  (1986)  but  was  modified  extensively  in  its  application  to  the 
Chesapeake  Bay  study  (Johnson  et  al.  1991).  These  modifications  include 
implementing  different  basic  numerical  formulations  of  the  governing  equa¬ 
tions  as  well  as  substantial  recoding  to  obtain  a  more  computationally  efficient 
model.  Physical  processes  impacting  circulation  which  can  be  modeled 
include  tide,  wind,  river  inflow,  and  the  effect  of  the  earth’s  rotation  (i.e., 
Coriolis  effect). 

A  key  attribute  of  CH3D  is  its  ability  to  define  a  basin  in  a  boundary-fitted 
coordinate  system,  allowing  grid  coordinate  lines  to  conform  with  irregular 
coastal  features,  such  as  a  shoreline  or  navigation  channel.  The  solution  algo¬ 
rithm  employs  an  external-internal  mode-splitting  technique.  In  the  external 
mode,  finite  difference  approximations  of  the  vertically-integrated  Navier- 
Stokes  equations  are  solved,  yielding  water  surface  elevations  and  depth- 
averaged  x-  and  y-directed  unit  flow  rates.  This  information  is  then  processed 
in  the  internal  mode  to  determine  the  x-,  y-,  and  z-directed  velocity  distribu¬ 
tions  through  the  water  column.  Because  the  Indian  River-Rehoboth  Bay 
model  application  is  a  two-dimensional,  depth-averaged  (i.e.,  external)  mode 
application,  the  internal  mode  is  not  discussed  in  this  chapter. 

This  chapter  describes  the  governing  equations  used  in  CH3D  and  the  for¬ 
mulation  of  bottom  and  surface  shear  stress  terms  as  well  as  the  Coriolis 
effect  The  process  of  transforming  and  non-dim ensionalizing  the  governing 
equations  is  discussed  and  finite -difference  approximations  of  the  governing 
equations  are  presented. 


Governing  Equations 

The  hydrodynamic  equations  used  in  CH3D  are  derived  from  the  classical 
Navier-Stokes  equations  formulated  in  a  Cartesian  coordinate  system 
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(Figure  4-1).  Assuming  that  the  vertical  water  accelerations  are  small  in  com¬ 
parison  with  the  gravitational  acceleration  (i.e.,  hydrostatic  pressure  conditions 
exist)  and  that  the  fluid  is  homogeneous  and  incompressible,  the  depth- 
averaged  approximation  yields  the  following  in-plan,  two-dimensional  form  of 
the  governing  equations: 

x-Momentum 
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Continuity 
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where 

x,y,t=  independent  space  and  time  variables 

U,  V  =  unit  flow  rate  components  in  the  x-  and  y -directions,  respectively 
H  -  total  water  depth  (h+S) 

S  =  water  surface  displacement  measured  relative  to  an  arbitrary 
datum 

h  =  static  water  depth  measured  from  the  same  datum 
g  =  gravitational  acceleration 
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Figure  4-1 .  Cartesian  coordinate  system  definition  sketch 


/=  Coriolis  parameter 

xsr  xsy  -  surface  shear  stress  in  the  x-  and  y -directions,  respectively 
iBx,  Xgy  *  bottom  shear  stress  in  the  x-  and  y -directions,  respectively 
p  =  water  density  (assumed  to  be  constant) 

Ah  =  generalized  dispersion  coefficient 


p  =  pressure 

Bottom  Shear  Stress  Formulation 


CH3D  uses  the  following  quadratic  expression  to  represent  the  bottom  shear 
stress  in  die  x-momentum  equation: 


xBx 


g  Jn 2  *  t/2  u 

c)h2 


(4) 


where 

Cz  =  Chezy’s  resistance  factor 
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A  similar  expression  is  used  for  xBy  in  the  y -momentum  equation. 

Rather  than  specifying  the  Chezy  resistance  factor,  Manning’s  n,  which  is 
independent  of  depth,  is  input  to  the  model.  These  coefficients  are  related 
through  the  following  equation: 


n 


(5) 


In  addition,  CH3D  has  an  option  for  designating  Manning’s  n  as  a  function  of 
depth  to  specify  changes  in  bottom  roughness  at  different  water  depths. 


Surface  Shear  Stress  Formulation 

The  surface  shear  stress  x,  is  formulated  as: 


t,  -  P.  CD\W\W 


(6) 


where  p«  is  the  a.r  density,  W  is  the  wind  velocity,  and  CD  is  a  dimensionless 
wind  drag  coefficient  CH3D  uses  the  wind  drag  formulation  presented  in 
Garratt  (1977): 

c  =  (0.75  +  0.067(0)  (7] 


where  ©  is  the  wind  speed.  This  formulation  requires  wind  speeds  specified  in 
units  of  meters  per  second.  An  upper  limit  of  3.0xlff3  is  applied  to  this  coeffi¬ 
cient  Thus,  for  wind  speeds  greater  than  63  knots,  a  constant  drag  coefficient 
is  applied. 


Coriolis  Effect 

Although  it  is  not  a  true  force,  foe  Coriolis  effect  accounts  for  the  apparent 
deflection  in  a  fluid’s  trajectory  that  is  induced  by  the  rotation  of  foe  Earth. 
The  Coriolis  parameter/ is  expressed  as: 


/  *  2vsinX 


(8) 
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where  v  is  the  angular  speed  of  the  Earth's  rotation  (7.292  x  10" 5  rad/sec)  and 
X  is  the  latitude  of  the  study  area  (38  deg  36.5  min). 

Non-Dimensionalization  of  Governing  Equations 

The  dimensionless  form  of  the  governing  equations  are  used  to  facilitate 
relative  magnitude  comparisons  of  the  various  terms  in  the  governing  equations 
and  to  minimize  the  effects  of  round-off  errors  during  computations.  The 
following  dimensionless  variables  are  used: 

(u\  v  \  w  *)  =  (u,  v,  wXyZ^AJ, 

(x\  y  *,  z  *)  =  (x,  y,  zXJZyX, 

(<•  <,)  =  (C  TyVp#fZUr 
t*  =  tf 

S*  «  gS/fUrXr  -  s/s, 

Ah  =  Ah^ht 

k;  « 

These  definitions  yield  the  following  dimensionless  parameters  in  the  govern¬ 
ing  equations: 

Froude  Number  F,  *  U/CgZ,)"2 

Rossby  Number  Ro  *  U/fX, 

Densimetric  Froude  Number  FrD  *  F,  /  Jz 

where 

e  =  (pr  -  pj/p„ 


Ur ,  p, ,  Xr ,  Z,  •  Ah,  ,  and  are  arbitrary  reference  values  of  the  velocity, 
density,  length,  depth,  dispersion,  and  diffusion. 

Using  the  dimensionless  variables  (asterisks  have  been  dropped)  and  the 
parameters  previously  defined,  the  vertically  integrated  equations  constituting 
the  external  mode  are: 
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Transformation  of  Governing  Equations 

The  governing  equations  contained  in  CH3D  were  developed  in  a  nun- 
orthogonal  curvilinear  coordinate  system.  This  system,  however,  necessitates 
the  transformation  of  the  governing  equations  into  a  non-orthogonal 
curvilinear  or  boundary-fitted  coordinate  system  (£,q).  Both  independent  (e.g., 
x,  y )  and  dependent  variables  (e.g.,  (/,  V)  in  the  governing  equations  are 
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transformed  into  the  (^,rj)  curvilinear  system.  CH3D  employs  contravariant 
components  (as  opposed  to  covariant  components)  in  the  transformation  of  the 
governing  equations,  therefore  velocities  are  defined  perpendicular  to  a  cell 
face,  as  opposed  to  parallel  to  a  cell  face. 

The  flow  rate  components  in  physical  space  (i.e.,  U{i)  and  V(  j))  are  related 
to  the  contravariant  components  (i.e.,  U\  V\  UK  Vj)  by  the  following 
equations: 


1/(0 


(12) 


VO) 


OULu*  +  J!?Lvj 
1ST  TiT 


(13) 


where  gy  is  the  metric  tensor  defined  as: 


*  +  *  Vn +  yy* 

v? +  v<  4 +  y\ 


(14) 


or 


&ll  &12 
#21  &22 


(15) 


and  |  g  |  is  the  determinant  of  the  metric  tensor  gy: 


I#  I  *  &11 8j2  ~ 


(16) 


Whereas  scalar  quantities  in  the  physical  plane  are  identical  in  the  transformed 
plane,  all  spatial  derivatives  containing  these  terms  must  be  transformed.  The 
surface  slope  terms  are  transformed  as  follows: 
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where  g M  are  inverse  metric  tensor  components: 
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or 


The  transfonned  governing  equations  developed  by  Sheng  (1986)  are  as 
follows.  The  inertial  and  diffusion  terms  in  contravariant  coordinates  are  quite 
lengthy  and  thus  are  omitted  in  this  report  However,  these  terms  are  pre¬ 
sented  in  Johnson  et  aL  (1991). 
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where 


U,V  =  contravariant  unit  flow  rate  components  in  die  transformed 
plane  (superscripts  have  been  dropped  for  convenience) 

I  gs  I  *  determinant  of  the  metric  tensor,  |g|  ,  at  an  S-point 
(Figure  4-2) 

I  ga  I ,  I  gv  I  =  determinant  of  the  metric  tensor,  |  g  |  ,  at  a  U-face  and  V-face, 
respectively 

U  =  average  x-direction  unit  flow  rate  at  a  V-face 
V  =  average  y -direction  unit  flow  rate  at  a  t/-face 

Finite  Difference  Approximations 
of  Governing  Equations 

The  finite  difference  approximations  to  the  governing  equations  are  based 
on  a  Eulerian  system  where  the  velocities  and  water  surface  fluctuations  are 
computed  at  discrete  locations  within  tire  flow  field.  A  network  of  grid  cells  is 
used  to  define  the  parameter  locations.  A  representative  grid  cell  in  compu¬ 
tational  space  (£,t0  is  shown  in  Figure  2.  In  this  staggered  grid,  the  water 
surface  fluctuation  is  defined  at  the  cell  center  (ij),  ^-direction  unit 
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flow  rates  (10  are  defined  at 
the  “west”  (ij)  and  “east” 
(i+l-j)  cell  faces,  and  the 
ri -direction  unit  flow  rates 
(VO  are  computed  at  the 
“south"  (ij)  and  “north” 
(ij+l)  cell  faces.  The  finite 
difference  approximations  of 
the  governing  equations 
follow.  Note  that  the  conti¬ 
nuity  equation  is  split  into 
two  parts.  The  sum  of  these 
equations  is  the  original 
continuity  equation. 

Figure  4-2.  Variable  positions 
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where 

n  =  previous  time  level 
*  =  intermediate  time  level 
n+1  =  solve  for  this  time  level 
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6  =  weighting  factor  between  successive  time  levels 
/  =  inertia 
FRC  =  friction 
D  -  dispersion 
n-Momentum 
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The  computational  procedure  used  in  CH3D  is  based  on  an  Alternating 
Direction  Implicit  scheme  (Roache  1976).  Using  this  method,  the  and 
n -momentum  equations  are  solved  separately,  and  each  calculation  in  time  is 
made  in  two  stages.  In  the  first  stage,  the  ^-continuity  and  ^-momentum  equa¬ 
tions  are  solved  along  each  row  of  the  grid  to  progress  from  time  level  n  to  an 
inteimediate  time  level  *,  indicated  by  superscript  *.  The  £ -direction  unit  flow 
rate  components  and  water  surface  fluctuations  are  solved  implicitly,  and  the 
T| -direction  unit  flow  rate  components  are  supplied  from  time  level  n.  The 
^-direction  unit  flow  rates  from  this  step  represent  those  at  time  level  n+1  , 
whereas  the  water  surface  fluctuations  are  only  an  approximation  to  those  at 
time  level  n+1.  The  r| -direction  unit  flow  rate  components  remain  at  time 
level  n.  In  the  second  stage,  the  q-continuity  and  q-momentum  equations  are 
solved  along  each  column  for  the  q-direction  unit  flow  rates  and  the  water 
surface  fluctuations  at  time  level  n+1.  ^-direction  unit  flow  rate  components 
are  supplied  from  the  first  stage  calculations. 

As  shown  in  the  finite  difference  approximations  to  the  governing  equations, 
a  weighting  factor  8  is  used  to  place  the  water  surface  slope  and  bottom  fric¬ 
tion  terms  between  time  levels  n  and  n+1.  When  the  weighting  factor  equals 
0.0,  these  terms  are  evaluated  at  the  previous  time  level  n  (explicit  treatment), 
whereas  when  the  weighting  factor  equals  1.0,  they  are  evaluated  at  the  new 
time  level  n+1  (implicit  treatment).  Usually  a  value  between  0.0  and  1.0  is 
used. 
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Chapter  V:  Hydrodynamic 
Model  Application 


In  this  chapter,  die  application  of  CH3D  to  the  Indian  River  Bay-Rehoboth 
Bay  system  is  described.  A  discussion  of  numerical  grid  development  and  the 
bathymetry  used  to  represent  the  study  area  is  given.  The  model  calibration 
and  validation  procedure  is  presented,  including  a  discussion  of  boundary 
forcing  conditions  used  in  the  model  and  an  analysis  of  model  results. 


Numerical  Grid  Development 

CH3D  uses  a  non-orthogonal,  boundary-fitted  grid  to  represent  the  limits  of 
a  study  area.  Using  a  boundary-fitted  grid  allows  coordinate  lines  to  conform 
to  the  irregular  boundaries  of  land  masses.  In  this  way,  the  delineation 
between  each  bay  and  the  surrounding  land  can  be  made.  Similarly,  the  delin¬ 
eation  between  an  island  and  the  surrounding  water  can  be  made.  Note  that 
CH3D  does  not  contain  an  adaptive  gridding  algorithm,  therefore  there  is  no 
means  to  allow  for  flooding  and  drying  of  low-lying  areas.  To  maintain  the 
appropriate  water  volume  in  die  system,  active  cell  volumes  in  low-lying  areas 
were  increased. 

The  boundary-conforming  grid  used  for  the  Indian  River  Bay-Rehoboth  Bay 
hydrodynamic  study  was  constructed  using  program  Eagle  (Thompson  1985). 
The  73  x  68  cell  grid  generated  by  Eagle  contains  approximately  2000  active 
cells  with  cells  concentrated  through  the  inlet  and  Middle  Island  areas  (Fig¬ 
ure  5-1).  Cell  sizes  through  the  inlet  were  approximately  30-m  by  120-m  in 
the  north-south  and  east-west  directions,  respectively.  Cell  sizes  near  the  Mid¬ 
dle  Island  area  were  as  small  as  60-m  by  120-m  in  the  east-west  and  north- 
south  directions,  respectively.  A  1988  hydrographic  survey  of  the  Indian  River 
Bay-Rehoboth  Bay  system  conducted  by  the  U.S.  Army  Engineer  District, 
Philadelphia  provided  the  necessary  bathymetric  information.  This  information 
was  translated  to  a  depth  in  each  grid  cell  by  applying  a  3-point  linear  interpo¬ 
lation  scheme  to  die  data.  Depths  in  the  dual  bay  system  arc  generally  less 
than  2  m  mean  low  water  (mlw). 
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Figure  5-1.  Indian  River-Rehoboth  Bay  grid 


Model  Calibration  and  Validation 


The  purpose  of  the  calibration  procedure  is  to  insure  that  the  model  accu¬ 
rately  predicts  hydrodynamic  conditions  within  a  given  study  area.  The  accu¬ 
racy  of  model  results  is  greatly  influenced  by  the  accuracy  of  boundary  forcing 
conditions,  representation  of  the  geometry  of  the  study  area  (i.e.  bathymetry 
and  land/water  interface),  and  to  a  lesser  degree,  the  choice  of  certain  “calibra¬ 
tion"  parameters  (Mark  et  al.  1993).  In  the  calibration  procedure,  these  para¬ 
meters  are  adjusted  to  obtain  the  best  agreement  between  model  results  and 
measured  field  data.  In  essence,  if  the  study  area  is  accurately  depicted  with 
the  numerical  grid  and  die  boundary  conditions  are  correct,  model  results  will 
be  “in  the  ballpark"  and  the  small  adjustments  to  calibration  parameters  are  a 
means  of  “fine-tuning”  model  results. 

The  purpose  of  the  validation  procedure  is  to  confirm  that  the  calibrated 
model  can  replicate  hydrodynamic  conditions  independent  of  the  conditions  for 
which  it  was  calibrated.  In  this  procedure,  the  model  is  applied  without  any 
adjustment  to  the  calibration  parameters.  A  good  comparison  between  model 
results  and  measured  data  in  the  validation  process  means  that  the  model 
should  accurately  simulate  hydrodynamics  for  other  time  periods. 

The  time  periods  selected  for  model  calibration  and  validation  for  the 
Indian  River  Bay-Rehoboth  Bay  system  were  time  periods  when  there  was 
available  field  data  for  comparison  to  model  results  with  sufficient  length  of 
record  and  sufficient  spatial  extent  to  cover  the  model  domain.  Another  factor 
influencing  the  time  periods  selected  was  the  fact  that  the  hydrodynamic  model 
results  were  used  to  supply  flow  conditions  to  the  water-quality  (WQ)  model. 
Therefore,  the  hydrodynamic  model  should  be  tested  over  a  time-span  consis¬ 
tent  with  the  transport  of  contaminants  within  this  system.  During  a  low  flow 
time  period,  river  discharges  may  lack  sufficient  momentum  to  flush  contami¬ 
nants  from  the  system.  During  extreme  events,  large  oscillations  in  the  water 
surface  and/or  high  river  discharges  can  provide  sufficient  momentum  to  trans¬ 
port  contaminants.  Because  of  the  importance  of  these  factors,  the  model  was 
calibrated  and  validated  for  a  cycle  of  7  days. 

Calibration  Overview 

CH3D  was  calibrated  for  the  Indian  River  Bay-Rehoboth  Bay  study  area  for 
the  one-week  time  period  29  June  -  5  July  1988.  Both  water  surface  elevation 
and  velocity  data  were  only  available  for  this  time  period.  Water  surface 
elevation  data  were  also  available  at  all  tide  gauges  from  1  October  - 
31  December  1988  and  comparisons  were  made  for  one  week  during  this  time 
period  (14  -  20  October  1988)  as  pan  of  the  validation  process. 

The  open  boundary  in  Indian  River  Inlet  was  driven  with  a  time  series  of 
water  surface  elevation  recorded  by  the  tide  gauge  adjacent  to  die  U.S.  Coast 
Guard  Station  (T5  in  Figure  2-1).  The  grid  boundary  location  was  selected  to 
coincide  with  the  gauge  location.  However,  grid  cells  in  this  region  are 
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approximately  120  m  wide  and  therefore  “precise”  placement  is  considered  to 
be  within  60  m  of  the  actual  gauge  location.  In  addition,  the  tide  gauge  is 
attached  to  a  pier  in  the  Coast  Guard  Boat  Basin  and  is  therefore  somewhat 
sheltered.  Because  the  influence  of  these  factors  is  not  known,  water  levels 
were  not  adjusted  in  phase  or  amplitude  to  account  for  these  factors.  As  a 
consequence,  water  surface  levels  and  velocities  computed  by  the  model  may 
contain  phase  errors  of  several  minutes. 

Tide  data  were  collected  at  the  USCG  Station  at  15 -minute  intervals,  how¬ 
ever,  water  levels  were  supplied  to  the  model  at  a  1-hr  interval  using  a  lowlow 
pass  (LLP)  filter  on  the  recorded  time  series.  Water  levels  were  lineariy- 
inteipolated  by  the  model  for  each  time  step  between  the  1-hr  input  values  and 
were  uniformly  assigned  to  every  cell  on  the  open  boundary.  No  water  surface 
level  gradient  was  imposed  along  this  boundary. 

Note  that  after  completing  the  initial  calibration  procedure,  a  simulation  was 
made  using  tire  original  15-min  interval  tidal  data  for  the  inlet  boundary  forc¬ 
ing  function  and  the  results  were  not  altered  significantly.  Therefore,  using  the 
1-hr  filtered  data  was  appropriate. 

The  time-history  of  wind  speed  and  directions  used  in  die  model  were 
recorded  by  the  Delaware  Power  and  Light  (DPL)  wind  anemometer.  These 
data  were  selected  over  wind  data  at  the  “old”  USCG  station  because  the  DPL 
record  was  more  complete.  Short  time  periods  missing  from  the  DPL  record 
were  linearly-interpolated  between  recorded  values.  DPL  data  for  30  -  31 
October  1988  were  extremely  high,  which  when  used  as  input  to  the  hydrody¬ 
namic  model  caused  grid  cells  to  become  dry,  thus  leading  to  model  instabili¬ 
ties.  After  consulting  with  DPL  personnel,  it  was  determined  that  these 
extremely  high  wind  velocities  were  inappropriate  and  were  adjusted  to  more 
reasonable  estimates  for  that  time  period.  Wind  data  were  supplied  to  the 
model  at  1-hr  intervals  and  were  linearly-interpolated  for  time  steps  between 
tire  1-hr  input  values. 

Discharge  data  were  available  from  1  January  -  30  September  1988  for 
Vines  Creek  and  Millsboro  Pond  (at  the  head  of  Indian  River).  Flows  at 
remaining  locations  throughout  the  system  (Figure  3-1)  were  estimated  by  the 
ratio  technique  described  in  Chapter  3.  Estimation  techniques  were  also  used 
to  complete  the  Vines  Creek  and  Millsboro  Pond  datasets  outside  die  period  of 
record.  Discharge  data  were  provided  by  USGS  as  daily  average  flow  rates, 
therefore  model  input  consisted  of  one  flow  rate  for  each  day  of  the  simula¬ 
tion.  Flow  rates  were  updated  at  each  time  step  of  the  simulation  by  linearly 
interpolating  between  daily  values. 

Data  available  for  comparison  with  model  results  included: 

1)  a  time  series  of  water  surface  elevations  recorded  at  each  of  the 
following  USGS  stations  (Figure  2-1): 
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a)  Vines  (Tl) 

b)  Fomets  (T2) 

c)  Massey  (T3) 

d)  Dewey  (T4) 

Water  surface  elevation  data  were  collected  at  15-minute  intervals.  Water 
level  data  collected  at  the  USCG  station  were  used  to  drive  the  inlet  boundary 
condition  and  therefore  could  not  be  used  for  comparison  to  model  results. 

2)  currents  speeds  without  directions  collected  at  four  locations  near 
Middle  Island  (CM1,  CM2,  CM3,  and  CM4  in  Figure  2-1). 


Calibration  Conditions 


The  hydrodynamic  model  was  used  to  simulate  the  year  1988  so  that  the 
WQ  team  could  begin  simulations  with  the  water  quality  model.  The  actual 
calibration  period,  however,  was  only  153  hours  (29  June  1500  -  5  July  2400). 
Data  were  available  at  all  four  tide  gauges  and  current  data  were  also  collected 
for  25  hrs  at  the  4  velocity  stations  (30  June  hr  0600  -  1  July  hr  0700)  during 
this  time  period.  Data  at  die  four  tide  gauges  showed  a  variation  in  range 
from  a  low  of  0.6  m  at  the  Dewey  gauge  to  a  high  of  1.1  m  at  the  Vines 
Creek  gauge.  River  discharge  for  the  calibration  time  period  was  extremely 
low.  The  Millsboro  Dam  flow  rate  was  0.76-0.91  m3/sec  and  the  Vines  Creek 
flow  rate  was  less  than  0.03  m3/sec.  Peak  flow  rates  for  1988  were  approxi¬ 
mately  5.66  m3/sec  at  Millsboro  Dam  and  2.83-4.25  m3/sec  at  Vines  Creek. 
Winds  during  the  calibration  period  were  1.3-5. 8  m/sec  generally  from  a  north¬ 
erly  direction. 


Calibration  Procedure 

Calibration  was  achieved  primarily  through  adjustments  to  the  bottom  fric¬ 
tion  coefficient.  Manning’s  n.  Initially  a  global  bottom  friction  coefficient  was 
specified  throughout  the  grid  and  the  model  was  run  repeatedly  with  different 
global  values  (0.020-0.035)  to  determine  an  approximate  value  of  n.  The  best 
global  value  was  determined  by  comparing  model  predicted  water  surface 
elevations  to  measured  water  surface  elevations  at  the  four  tide  stations 
described  previously.  To  achieve  the  more  difficult  task  of  replicating  water 
velocities,  the  global  bottom  friction  coefficient  was  refined  and  tested  over  a 
narrower  range  of  values.  In  the  last  stage,  Manning’s  n  values  were  adjusted 
locally  in  shallow  regions  to  account  for  increased  friction  drag  in  these  areas. 
The  best  comparison  between  model  results  and  measured  data  was  achieved 
with  a  global  Manning’s  n  of  0.025.  Shallower  areas  (less  than  .3  m)  were 
assigned  a  friction  value  of  0.045. 
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Analysis  of  Calibration  Results 

Water  Surface  Elevation 

The  model  accurately  reproduced  the  amplitude  of  the  water  surface  level 
time  histories  recorded  at  the  Vines  Creek,  Potnets,  Massey,  and  Dewey 
gauges,  however,  the  model  consistently  leads  the  prototype  data  by  30  min  at 
all  gauges  (Figures  3-2  through  5-5).  The  modeled  tide  travels  too  quickly 
from  the  inlet  to  die  gauges.  Scatter  plots  of  predicted  versus  measured  water 
surface  levels  show  the  typical  elliptic  shape  indicative  of  a  phase  shift  (Fig¬ 
ures  5-6  through  5-9).  The  root-mean-square  (rms)  error  at  the  Vines,  Potnets, 
Massey,  and  Dewey  gauges  was  computed  to  be  7.1,  8.3,  3.0,  and  2.7  cm, 
respectively.  A  thorough  search  for  the  cause  of  this  discrepancy  led  to  many 
interesting  possibilities. 

First,  water  surface  level  data  computed  by  the  model  is  written  to  an  out¬ 
put  file  at  an  interval  of  30  minutes  -  the  exact  phase  discrepancy.  Although 
this  seemed  to  be  a  likely  source  for  the  discrepancy,  it  was  not  die  cause  in 
this  case.  The  model  was  run  with  an  output  interval  of  15  minutes,  and  the 
results  were  identical  to  the  previous  results.  Another  possibility  for  the  dis¬ 
crepancy  was  data  interpretation,  i.e.,  a  plotting  error.  Rather  than  plotting  the 
results,  the  actual  values  produced  by  die  model  were  evaluated  and  compared 
with  the  prototype  values,  and  the  discrepancy  was  real. 

Interestingly,  all  four  gauges  are  consistently  out  of  phase  with  the  proto¬ 
type,  regardless  of  their  distance  from  tbs  inlet.  Therefore,  the  error  lies  either 
with  the  boundary  condition  imposed  at  the  inlet  or  between  the  inlet  boundary 
and  die  closest  gauge  (Massey  (T3)).  First  the  inlet  boundary  forcing  function 
was  investigated  as  the  source  for  the  discrepancy.  The  grid  boundary  location 
was  specifically  selected  to  coincide  with  die  USCG  station  gauge  location 
(T5).  However,  grid  cells  in  this  region  are  approximately  120  m  wide  and 
therefore  ‘‘precise”  placement  is  considered  to  be  within  60  m  of  the  actual 
gauge  location,  in  addition,  the  tide  gauge  is  attached  to  a  pier  in  the  Coast 
Guard  Boat  Basin  and  is  therefore  somewhat  sheltered.  The  influence  of  these 
factors  may  contribute  to  the  phase  discrepancy  observed  in  the  water  surface 
levels  computed  by  the  model,  however,  die  contribution  would  be  minor. 

Tide  data  were  collected  at  the  USCG  Station  using  a  digital  recorder  with 
a  15 -minute  punch  interval,  however,  water  levels  were  supplied  to  the  model 
at  a  1-hr  interval  using  a  low-low  pass  (LLP)  filter  on  die  recorded  time  series. 
Water  levels  were  linearly-interpolated  by  die  model  for  each  time  step 
between  the  1-hr  input  values  and  were  uniformly  assigned  to  every  cell  on  the 
open  boundary.  Could  die  1-hr  interval  between  boundary  forcing  data  be  the 
cause  for  the  phase  lag?  To  investigate  this,  die  model  was  rerun  with  an  inlet 
boundary  forcing  interval  of  15  minutes.  The  results  were  not  altered  signifi¬ 
cantly,  therefore  die  1-hour  forcing  interval  was  appropriate. 
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Figure  5-2.  Tidal  calibration  at  Vines  Creek  gauge 


Figure  5-3.  Tidal  calibration  at  Pot  nets  gauge 
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Figure  5-5.  Tidal  caltoraiion  at  Dewey  Beach  gauge 
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Figure  5-6.  Scatter  plot  of  tidal  calibration  at  Vines  Creek  gauge 


Figure  5-7.  Scatter  plot  of  tidal  calibration  at  Potnefs  gauge 
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Figure  5-9.  Scatter  plot  of  tidal  caHbration  at  Dewey  Beach  gauge 
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The  choice  of  a  water  surface  level  boundary  condition  was  examined. 
Choosing  to  force  the  model  at  the  inlet  boundary  with  water  surface  elevation 
data  was  dictated  by  the  considerable  quantity  of  water  surface  level  data  avai¬ 
lable  there.  However,  velocities  through  Indian  River  Inlet  are  significant  (1.5- 
2.5  m/sec).  If  the  model  was  forced  with  flow  through  the  open  boundary 
rather  than  a  water  level  time  series,  the  boundary  condition  would  include  the 
velocity  head  contribution.  The  relative  contribution  of  elevation  head  and 
velocity  head  may  not  be  computed  correctly  by  merely  specifying  a 
water  surface  elevation  boundary  condition.  This  factor  could  have  an  impact 
on  phasing  of  the  tides.  Therefore,  it  is  concluded  that  the  location  of  the  inlet 
boundary  precludes  the  model  from  simulating  the  large  head  loss  through  the 
inlet.  The  boundary  condition  applied  across  the  inlet  boundary  does  not 
include  the  velocity  head  contribution  nor  its  interaction  with  the  elevation 
head.  These  factors  most  probably  caused  the  observed  phasing  discrepancy. 

Regardless  of  the  cause  of  the  phase  shift,  it  is  insignificant  from  a  water 
quality  standpoint.  By  adjusting  the  boundary  forcing  function  30  minutes,  the 
model  results  are  in  phase  with  the  prototype  data  at  all  four  gauges  (Fig¬ 
ures  5-10  through  5-13).  In  addition,  the  scatter  plots  of  predicted  versus 
measured  water  surface  level  lose  their  elliptic  shape  and  show  a  linear  rela¬ 
tionship  between  model  results  and  prototype  data  (Figures  5-14  through  5-17). 
The  ims  error  at  the  Vines,  Potnets,  Massey,  and  Dewey  gauges  was  reduced 
to  3.1,  4.3,  2.0,  and  2.2  cm,  respectively. 


Velocity 

Only  25  hours  of  velocity  data  were  available  for  comparison  with  model 
results.  All  four  velocity  stations  (CM1  through  CM4)  are  in  the  vicinity  of 
Middle  Island  (Figure  2-1).  The  model  overpredicts  the  water  velocity  east  of 
Middle  Island  (at  CM1  and  CM2),  is  fairly  accurate  west  of  Middle  Island  at 
CM3,  and  underpredicts  flood  velocities  and  overpredicts  ebb  velocities  at 
CM4  (Figures  5-18  through  5-21). 


Validation  Conditions 

CH3D  was  validated  for  the  Indian  River  Bay-Rehoboth  Bay  study  area  for 
the  one-week  time  period  14  -  20  October  1988.  Water  surface  elevation  data 
were  available  at  the  same  four  tide  gauges  referred  to  in  the  calibration,  for 
the  validation  time  period.  A  comparison  of  model  results  to  these  data  was 
used  to  validate  CH3D.  Data  at  the  four  tide  gauges  showed  a  variation  in 
range  from  a  low  of  0.6  m  at  Dewey  gauge  and  a  high  of  1.1  m  at  the  Vines 
gauge. 

As  in  the  calibration  process,  the  open  boundary  at  Indian  River  Inlet  was 
driven  with  a  time  series  of  water  surface  elevation  data  recorded  by  the  tide 
gauge  adjacent  to  the  USCG  Station  (T5  in  Figure  2-1).  Tide  data  were  col¬ 
lected  at  the  USCG  Station  at  15-minute  intervals,  however,  water  levels  were 
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Figure  5-10.  Tidal  calibration  at  Vines  Creek  gauge 


Figure  5-11.  Tidal  caHbration  at  Potnets  gauge 
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Figure  5-12.  Tidal  calibration  at  Massey’s  Ditch  gauge 


Figure  5-13  Tidal  calibration  at  Dewey  Beach  gauge 
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Figure  5-14.  Scatter  plot  of  tidal  calibration  at  Vines  Creek  gauge 
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Figure  5-15.  Scatter  piot  of  tidal  calibration  at  Potnets  gauge 
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Figure  5-16.  ScaRer  plot  of  tidal  calibration  at  Massey’s  Ditch  gauge 
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Figure  5-17.  Scatter  plot  of  tidal  caUbratlon  at  Dewey  Beach  gauge 
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Figure  5-18.  Velocity  calibration  at  Middle  Island  Station  1 


Figure  5-19.  Velocity  calibration  at  Middle  Island  Station  2 
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Figure  5-20.  Velocity  calibration  at  Middle  Island  Station  3 


Figure  5-21.  Velocity  calibration  at  Middle  Island  Station  4 
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supplied  to  the  model  at  a  1 -hr  interval  using  a  low-low  pass  filter  on  the 
recorded  time  series.  Water  levels  were  linearly  interpolated  by  the  mode)  for 
each  time  step  between  the  1-hr  input  values  and  were  uniformly  assigned  to 
every  cell  on  the  open  boundary.  No  water  surface  level  gradient  was  imposed 
along  this  boundary. 

Wind  speed  and  direction  used  in  the  model  were  recorded  by  the  Delaware 
Power  and  Light  (DPL)  anemometer.  Wind  data  were  supplied  to  the  model  at 
1-hr  intervals  and  were  linearly  interpolated  for  time  steps  between  the  1-hr 
input  values. 

Discharge  data  for  the  validation  period  were  estimated  by  ratio  techniques 
described  in  Chapter  3.  River  discharge  rates  for  the  validation  period  were 
average.  The  Millsboro  Dam  flow  rate  was  approximately  1.66  m3/sec  and  the 
Vines  Creek  flow  rate  was  approximately  0.05  m3/sec.  Discharge  data  were 
estimated  as  daily  average  flow  rates,  therefore  model  input  consisted  of  one 
flow  rate  for  each  day  of  the  simulation.  Flow  rates  were  updated  at  each  time 
step  of  the  simulation  by  linearly  interpolating  between  daily  values. 


Analysis  of  Validation  Results 

As  with  the  calibration  period,  die  model  accurately  reproduces  the  ampli¬ 
tude  of  the  water  surface  level  time  histories  recorded  at  the  Vines  Creek, 
Potnets,  Massey,  and  Dewey  gauges,  however,  the  model  consistently  leads  the 
prototype  by  30  minutes  at  all  gauges  (Figures  5-22  through  5-25).  Scatter 
plots  of  predicted  versus  measured  water  surface  levels  show  the  typical  ellip¬ 
tic  shape  indicative  of  a  phase  shift  (Figures  5-26  through  5-29).  The  rms 
error  at  the  Vines,  Potnets,  Massey,  and  Dewey  gauges  was  computed  to  be 
5.1, 4.3,  3.3,  and  1.8  cm  respectively.  Again,  the  cause  of  the  discrepancy  is 
most  likely  due  to  die  boundary  condition  applied  across  the  inlet  boundary. 

By  adjusting  the  boundary  forcing  function  30  minutes,  die  model  results 
ate  in  phase  with  the  prototype  data  at  all  four  gauges  (Figures  5-30  through  5- 
33).  In  addition,  die  scatter  plots  of  predicted  versus  measured  water  surface 
level  lose  their  elliptic  shape  and  show  a  linear  relationship  between  model 
results  and  prototype  data  (Figures  5-34  through  5-37).  The  rms  error  at  die 
Vines,  Potnets,  and  Massey  gauges  was  reduced  to  2.8,  2.0,  and  1.3  cm, 
respectively.  The  rms  error  at  the  Dewey  gauge  did  not  change. 


Salinity  Overview 

The  salinity  distribution  in  the  Indian  River  Bay-Rehoboth  Bay  system 
varies  both  spatially  and  temporally.  Indian  River  Bay  has  a  strong  salinity 
gradient  in  the  east-west  direction,  whereas  Rehoboth  Bay  responds  fairly 
uniformly  (Figure  5-38).  The  strong  east-west  salinity  gradient  in  Indian  River 
Bay  is  clearly  the  result  of  saltwater  inflow  through  Indian  River  Inlet  (to  the 
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Figure  5-22.  Tidal  verification  at  Vines  Creek  gauge 


Figure  5*23.  Tidal  verification  at  Potnets  gauge 
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Figure  5-24.  Tidal  verification  at  Massey's  Ditch  gauge 
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Figure  5-25.  Tidal  verification  at  Dewey  Beach  gauge 
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Figure  5-26.  Scatter  plot  of  tidal  verification  A  Vines  Creek  gauge 


Figure  5-27.  Scatter  plot  of  tidal  verification  at  Pot  nets  gauge 
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Figure  5-29.  Scatter  plot  of  tidal  calixation  at  Dewey  Beach  gauge 
5-22 
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east)  and  freshwater  inflow  from  Indian  River  at  Millsboro  Dam  (to  the  west) 
(Kaipas  1978).  For  low  flow  conditions,  freshwater  does  not  penetrate  as  far 
into  Indian  River  as  under  high  flow  conditions,  resulting  in  a  lower  salinity 
gradient  from  the  inlet  boundary  to  the  river  boundary.  Conversely,  high  flow 
conditions  push  freshwater  far  into  Indian  River,  resulting  in  a  higher  salinity 
gradient  from  the  inlet  boundary  to  the  river  boundary.  Comparing  high  flow 
conditions  during  flood  and  ebb  shows  the  limits  of  salinity  excursion  in 
Indian  River  Bay  are  strongly  influenced  by  the  tide  (time  of  day).  For  high 
river  flow  conditions  and  flood  flow,  salinity  in  Indian  River  is  compressed  by 
tire  strong  river  flow  opposing  the  tidal  flood  flow  resulting  in  a  steep  salinity 
gradient  in  Indian  River  Bay.  Under  ebb  conditions,  the  salinity  gradient  is 
still  large,  but  the  excursion  limits  are  pushed  further  east  by  the  strong  river 
flow  and  ebbing  tide. 

Rehoboth  Bay  has  no  major  source  of  freshwater.  Because  of  this,  Reho- 
both  Bay  has  a  higher  salinity  concentration  than  most  of  Indian  River  Bay, 
and  Rehoboth  Bay  responds  nearly  uniformly  to  the  tidal  input  of  salinity. 
However,  there  is  a  weak  east-west  gradient  as  a  result  of  freshwater  inflow 
from  two  major  creeks  (Love  and  Herring)  on  the  western  shore  of  Rehoboth 
Bay  (Kaipas  1978)  (Figure  5-38). 

Temporal  variability  of  salinity  is  a  function  of  tidal  oscillation  (salt-water 
inflow)  and  variation  in  freshwater  inflow,  and  the  interaction  of  these  two 
factors.  Figures  5-39  and  5-40  show  the  time  variation  of  salinity  from  day 
260  to  day  265  of  1988  (September  16-20)  in  Indian  River  and  Indian  River 
Bay  (grid  cells  (9,29)  and  (30,24),  respectively),  and  Figures  5-41  and  5-42 
show  the  time  variation  of  salinity  near  the  center  of  Rehoboth  Bay  (grid  cells 
(42,54)  and  (44,60)  respectively). 

The  time  series  of  salinity  in  the  upper  reach  of  Indian  River  (Figure  5-39) 
shows  a  large  variation  in  salinity  with  tidal  oscillations.  Salinity  in  this 
region  ranges  from  a  low  of  13.7  ppt  to  a  high  a  21.5  ppt  with  an  average 
value  of  approximately  17-18  ppt.  The  important  point  is  that  salinity  varies 
nearly  8  ppt  during  one  tidal  cycle  in  the  upper  reach  of  Indian  River.  This  is 
critical  from  a  modeling  standpoint  because  of  the  sparseness  of  salinity  data 
for  comparison  to  model  results  and  the  missing  component  of  “time  of  day” 
in  the  salinity  data  collection  process. 

Figure  5-40  shows  that  the  time  variation  of  salinity  in  Indian  River  Bay 
has  less  influence  from  Indian  River  than  the  previous  location  and  is  therefore 
more  tidally  influenced.  Salinity  in  this  region  ranges  from  29.0  to  31.5  ppt 
with  an  average  value  of  approximately  30  ppt.  At  a  distance  from  the  influ¬ 
ence  of  freshwater  inflow,  the  variation  in  salinity  is  only  2.5  ppt. 

Figure  5-41  shows  that  salinity  in  the  south-central  portion  of  Rehoboth 
Bay  varies  slightly  with  the  tide,  having  a  range  of  only  1  ppt.  Salinity  in  this 
portion  of  Rehoboth  Bay  varies  between  30.5  and  31.5  ppt.  Figure  5-42  shows 
that  salinity  in  the  north-central  portion  of  Rehoboth  Bay  is  nearly  constant 
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Figure  5-30.  Tidal  verification  at  Vines  Creek  gauge 


Figure  5-31.  Tidal  verification  at  Potnets  gauge 
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Figure  5-32.  Tidal  verification  at  Massey’s  Ditch  gauge 


Figure  5-33.  Tidal  verification  at  Dewey  Beach  gauge 
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Figure  5-35.  Scatter  plot  of  tidal  verification  at  Pot  nets  gauge 
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Figure  5-36.  Scatter  plot  of  tidal  verification  at  Massey's  Ditch  gauge 
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Figure  5-37.  Scatter  plot  of  tidal  calibration  at  Dewey  Beach  gauge 
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Figure  5-38.  Indian  River  Bay/Rehoboth  Bay  salinity  distribution 


with  time  (31.2  ppt).  Rehoboth  Bay  has  no  major  freshwater  inflows  and  is 
therefore  tidaily-dom inated . 

In  conclusion,  Indian  River  Bay  is  the  more  complex  of  the  two  bays  in 
this  system  from  a  salinity  standpoint  Salinity  varies  both  spatially  and  tem¬ 
porally  in  Indian  River  Bay  due  to  the  interaction  of  salt  water  from  the  inlet 
and  freshwater  from  Indian  River  at  Mills!  o  Dam.  Rehoboth  Bay  is  nearly 
homogeneous.  There  is  no  major  source  oi  treshwater  inflow  to  Rehoboth  Bay 
and  it  is  therefore  clearly  tidally-dominated. 


Salinity  Calibration  and  Validation 

A  long-term  (1  year)  simulation  was  used  in  the  salinity  calibration.  The 
model  simulation  period  was  year  1988,  however  prototype  data  available  for 
this  time  period  was  limited  to  3  days  of  spatially-distributed  salinity.  A  time- 
series  of  salinity  at  one  location  would  therefore  be  limited  to  3  values  for  tire 
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Figure  5-39.  Temporal  variation  of  salinity  in  Indian  River 


Figure  5-40.  Temporal  variation  of  salinity  in  Indian  River  Bay 
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Figure  5-41.  Temporal  variation  of  salinity  in  Rehoboth  Bay  (42,54) 


Figure  5-42.  Temporal  variation  of  salinity  in  Rehoboth  Bay  (44,60) 
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entire  year.  Thus,  a  spatial  comparison  of  model  results  to  prototype  data  was 
used.  These  data  provide  a  picture  of  the  spatial  distribution  of  salinity  at 
14  stations  for  June  20,  July  25  and  September  19,  1988.  As  previously  dis¬ 
cussed,  it  is  important  to  note  that  the  exact  time  of  day  for  data  collection  is 
not  available.  This  factor  is  most  critical  in  Indian  River  where  the  temporal 
variation  of  salinity  can  be  as  much  as  8  ppt  during  one  tidal  cycle. 

For  illustration,  Figure  5-43  is  a  model  predicted  salinity  snapshot  of  the 
entire  study  area  at  12  noon  on  September  19,  1988  with  the  shaded  circles 
indicating  measured  salinity  values.  Note  that  measured  salinity  values  are 
single,  “instantaneous”  values,  however  the  time  of  day  that  the  data  were 
collected  is  not  known.  For  this  reason,  model  results  were  arbitrarily  saved  at 
12  noon.  (This  is  critical  in  areas  where  the  temporal  variability  of  salinity  is 
significant.)  Results  in  Rehoboth  Bay  are  excellent  Results  in  Indian  River 
Bay  are  very  good,  with  the  exception  of  the  upper  reach  of  Indian  River 
where  the  model  is  underpredicting  salinity.  The  model  salinity  saved  for  this 
region  is  in  the  range  14-17  ppt  however,  the  prototype  value  is  in  the  17- 
20  ppt  range. 

Reviewing  the  time-series’  of  salinity  for  the  central  portion  of  Rehoboth 
Bay  (Figure  5-41  and  5-42)  shows  that  salinity  is  nearly  constant,  therefore  the 
time  of  day  for  data  collection  is  not  critical.  In  contrast,  the  time-series  of 
salinity  in  the  upper  reach  of  Indian  River  (Figure  5-39)  shows  great  variability 
in  salinity  with  tide  (time  of  day).  Thus,  the  difference  between  model  results 
and  prototype  data  in  Indian  River  is  well  within  the  range  of  temporal  vari¬ 
ability  of  salinity  for  this  region. 

Although  the  model  results  are  within  the  range  of  temporal  variability, 
salinity  values  in  Indian  River  are  consistently  lower  than  the  prototype  values 
for  all  three  days  of  data  collection  in  1988.  Therefore,  several  techniques 
were  examined  to  obtain  a  better  salinity  calibration  in  the  upper  reach  of 
Indian  River.  The  technique  selected  was  to  increase  the  size  of  a  few  cells  in 
the  upper  reach  of  Indian  River.  The  added  cell  volume  decreased  the  velocity 
of  freshwater  inflow  thereby  allowing  salt  to  penetrate  further  into  Indian 
River.  In  addition,  a  spatially-variable  dispersion  coefficient  was  used  to  cali¬ 
brate  the  model.  A  value  of  12  m2/sec  was  used  in  Rehoboth  Bay  and  most  of 
Indian  River  Bay  with  the  exception  of  a  portion  of  the  middle  reach  of  Indian 
River  (cells  4-14).  A  spatially  variable  longitudinal  dispersion  coefficient  was 
used  to  calibrate  the  model  to  the  high  salinity  gradient  observed  in  this  region 
in  1988,  1989,  and  1990.  This  section  of  the  Indian  River  is  a  shallow  chan¬ 
nel  where  freshwater  from  the  river  meets  salt  water  from  the  lower  estuary  in 
addition  to  heated  water  discharged  from  the  DPL  power  plant  Smith  (1976) 
reported  that  in  the  presence  of  buoyant  concentration  either  from  a  river  or 
from  heated  water,  the  transverse  circulation  can  lead  to  a  marked  reduction  in 
longitudinal  dispersion.  It  is  believed  that  this  occurred  in  the  central  portion 
of  Indian  River.  For  this  reason,  a  smaller  dispersion  value  of  0.6  m2/sec  was 
used  in  cells  7  through  1 1  and  was  gradually  increased  to  12  m2/sec  in  both 
the  upstream  and  downstream  directions.  The  desired  salinity  distribution  was 
thus  obtained. 
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Figure  5-43.  Salinity  calibration  for  September  19, 1988 


A  second  long-term  (1-year)  simulation  was  used  in  the  salinity  validation. 
The  model  simulation  period  was  year  1989,  however  prototype  data  available 
for  this  time  period  was  limited  to  2  days  of  spatially-distributed  salinity.  A 
time-series  of  salinity  at  one  location  would  therefore  be  limited  to  2  values 
for  the  entire  year.  Thus,  a  spatial  comparison  of  model  results  to  prototype 
data  was  again  used.  These  data  provide  a  picture  of  the  spatial  distribution  of 
salinity  at  14  stations  for  April  24  and  October  17,  1989.  Figure  5-44,  depict¬ 
ing  the  April  24,  1989  model-prototype  comparison,  is  during  a  period  of  high 
flow  from  Indian  River.  The  results  are  excellent  with  the  exception  of  one 
location  in  Rehobo th  Bay.  At  gage  5,  the  prototype  data  indicates  a  higher 
salinity  (>29  ppt)  than  the  model  predicts  (26-29  ppt).  However,  the  actual 
salinity  value  is  only  29.9  ppt  and  the  discrepancy  is  therefore  small.  The 
lower  value  of  model  salinity  indicates  that  the  salinity  wedge  should  be  forced 
further  into  Rehoboth  Bay.  Figure  5-45,  depicting  the  October  17, 1989 
model-prototype  comparison,  is  during  a  period  of  low  flow  from  Indian  River. 
Model  results  are  excellent  at  all  gages. 
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Figure  5-44.  Salinity  validation  for  April  24, 1989 


Summer  Average  Salinity 

The  link  between  hydrodynamic  and  water  quality  modeling  is  in  the  con¬ 
servative  property  of  salinity.  In  water  quality  modeling,  summer  is  the 
chemically-  and  biologically-active  time  period.  Therefore,  a  comparison  of 
data  and  model  results  in  the  summer  months  is  important  to  ensure  the  accu¬ 
racy  of  concentration  transport  Since  the  data  available  in  the  summer  months 
precludes  a  time  series  comparison,  the  average  summer  conditions  were 
modeled  and  compared  to  all  measured  data  for  the  summers  of  1988,  1989, 
and  1990.  Summer  was  defined  in  the  model  as  July,  August,  and  September. 
Figures  5-46  through  5-48  show  comparisons  between  model  and  prototype 
data  for  the  summers  of  1988,  1989,  and  1990,  respectively.  The  three-month 
average  of  model  salinity  from  the  head  of  the  Indian  River  (Millsboro  Pond) 
to  Indian  River  Inlet  is  given  as  a  solid  line  and  prototype  data  at  8  salinity 
stations  on  three  specific  days  each  year  (July  24,  1989,  August  28,  1989, 
September  25,  1989;  July  16,  1990,  August  13,  1990,  September  24,  1990)  are 
given  as  +’s.  In  1988  only  two  days  of  salinity  data  (July  25,  September  19) 
were  available.  The  results  confirm  that  the  summer  trend  for  salinity  is 
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Figure  5-45.  Salinity  validation  for  October  17, 1989 


captured  by  the  model.  Some  discrepancy  between  model  and  data  is  evident 
near  Millsboro  Pond  in  1988  but  the  discrepancy  is  more  likely  due  to  sparse 
data  than  a  model  deficiency. 
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Chapter  VI:  Water  Quality 
Model  Kinetics 


Introduction 

The  water  quality  model,  dubbed  CE-QUAL-ICM,  originated  as  pan  of  a 
eutrophication  model  package  applied  to  Chesapeake  Bay  (Cerco  and  Cole 
1993).  Kinetic  formulations  from  the  Chesapeake  Bay  model  were  applied  to 
Indian  River-Rehoboth  Bay  with  minimal  modification.  Several  state  variables 
present  in  the  original  model  were  dropped,  however,  consistent  with  the  avail¬ 
able  data  base.  Some  simplifications  were  also  made,  largely  through  specifi¬ 
cation  of  parameter  values.  The  original  code  is  intact  and  additional  state 
variables  and  options  can  be  activated  in  die  future.  Active  state  variables  in 
Indian  River-Rehoboth  Bay  are  presented  in  Table  6-1.  The  remainder  of  the 
chapter  describes  the  state  variables,  details  kinetics  formulations,  and  lists 
coefficient  values. 


Table  6-1 

Water  Quality  Model  State  Variables 

Temperature 

Salinity 

Algae 

Dissolved  Organic  Carbon 

Labile  Particulate  Organic  Carbon 

Refractory  Paniculate  Organic  Cartoon 

Ammonium 

Nitrate 

Dissolved  Organic  Nitrogen 

Labile  Particulate  Organic  Nitrogen 

Refractory  Particulate  Organic  Nitrogen 

Total  Phosphate 

Dissolved  Organic  Phosphorus 

Labile  Particulate  Organic  Phosphorus 

Refractory  Particulate  Organic  Phosphorus 

Dissolved  Oxygen 

Algae 

Algae  are  a  primary  focus  of  the  eutrophication  model.  Algae  play  a  cen¬ 
tral  role  in  the  carbon,  nitrogen  and  phosphorus  cycles.  Dissolved  oxygen  is 
influenced  by  algal  production  and  respiration  in  the  water  column.  Light 
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extinction  is  paitiaily  dependent  attenuation  caused  by  algal  chlorophyll.  One 
algal  group  is  currently  activated  in  the  model. 


Organic  Cartoon 

Three  organic  carbon  state  variables  are  considered:  dissolved,  labile  partic¬ 
ulate,  and  refractory  particulate.  Labile  and  refractory  distinctions  are  based 
upon  the  time  scale  of  decomposition.  Labile  organic  carbon  decomposes  on  a 
time  scale  of  days  to  weeks  while  refractory  organic  carbon  requires  more 
time.  Labile  organic  carbon  decomposes  rapidly  in  the  water  column  or  the 
sediments.  Refractory  organic  carbon  decomposes  slowly,  primarily  in  the 
sediments,  and  may  contribute  to  sediment  oxygen  demand  years  after 
deposition. 


Nitrogen 

Nitrogen  is  first  divided  into  organic  and  mineral  fractions.  Organic  nitro¬ 
gen  state  variables  are:  dissolved  organic  nitrogen,  labile  particulate  organic 
nitrogen,  and  refractory  particulate  organic  nitrogen.  Two  mineral  nitrogen 
forms  are  considered:  ammonium  and  nitrate.  Both  are  utilized  to  fulfill  algal 
nutrient  requirements  although  ammonium  is  preferred  from  thermodynamic 
considerations.  The  primary  reason  for  distinguishing  the  two  is  that  ammon¬ 
ium  is  oxidized  by  nitrifying  bacteria  into  nitrate.  This  oxidation  can  be  a 
significant  sink  of  oxygen  in  the  water  column  and  sediments.  An  intermedi¬ 
ate  in  tire  complete  oxidation  of  ammonium,  nitrite,  also  exists.  Nitrite  con¬ 
centrations  are  usually  much  less  than  nitrate  and  for  modeling  purposes  nitrite 
is  combined  with  nitrate.  Hence  the  nitrate  state  variable  actually  represents 
the  sum  of  nitrate  plus  nitrite. 


Phosphorus 

As  with  carbon  and  nitrogen,  organic  phosphorus  is  considered  in  three 
states:  dissolved,  labile  particulate,  and  refractory  particulate.  Only  a  single 
mineral  form,  total  phosphate,  is  considered.  Total  phosphate  exists  in  two 
states  within  die  model  ecosystem:  dissolved  phosphate  and  phosphate  incorpo¬ 
rated  in  algal  cells.  Equilibrium  partition  coefficients  distribute  the  total 
among  two  states. 


Dissolved  Oxygen 

Dissolved  oxygen  is  required  for  the  existence  of  higher  life  forms.  Oxygen 
availability  determines  the  distribution  of  organisms  and  die  flows  of  energy 
and  nutrients  in  an  ecosystem.  Dissolved  oxygen  is  a  central  component  of  the 
water  quality  model. 
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Salinity 


Salinity  is  a  conservative  tracer  that  provides  verification  of  the  transport 
component  of  the  model  and  facilitates  examination  of  conservation  of  mass. 
Salinity  also  influences  the  dissolved  oxygen  saturation  concentration. 


Temperature 

Temperature  is  a  primary  determinant  of  the  rate  of  biochemical  reactions. 
Reaction  rates  increase  as  a  function  of  temperature  although  extreme  tempera¬ 
tures  result  in  the  mortality  of  organisms. 


Conservation  Of  Mass  Equation 

The  foundation  of  CE-QUAL-ICM  is  the  solution  of  the  mass-conservation 
equation  for  a  control  volume.  Control  volumes  in  the  water  quality  model 
correspond  to  one  or  more  cells  on  the  hydrodynamic  model  grid.  For  the 
Indian  River-Rehobo th  Bay  application,  a  two-dimensional,  vertically  inte¬ 
grated  form  of  the  mass-conservation  equation  is  employed.  CE-QUAL-ICM 
solves,  for  each  volume  and  for  each  state  variable,  the  conservation  of  mass 
equation: 


8V:  C;  " 


p-E<WE^*E 

j=l  j=l  °*j 


Si 


(6-1) 


=  volume  of  ith  control  volume  (m3) 

Cj  =  concentration  in  ith  control  volume  (gm  m'3) 

Qj  =  volumetric  flow  across  flow  face  j  of  ith  control  volume  (m3  sec1) 

Cj*  =  concentration  in  flow  across  flow  face  j  (gm  m*3) 

Aj  =  area  of  flow  face  j  (m2) 

Dj  =  diffusion  coefficient  at  flow  face  j  (m2  sec*1) 

n  =  number  of  flow  faces  attached  to  ith  control  volume 

St  =  external  loads  and  kinetic  sources  and  sinks  in  ith  control  volume 
(gm  sec*1) 

t,  x  =  temporal  and  spatial  coordinates 
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Solution  to  the  mass-conservation  equation  is  via  the  finite -difference 
method  using  the  QUICKEST  algorithm  (Leonard  1979). 

The  remainder  of  this  chapter  details  the  kinetics  portion  of  the  mass- 
conservation  equation  for  each  state  variable.  Terms  are  defined  where  they 
first  appear.  Coefficient  values  are  presented  in  a  table  at  the  end  of  the  chap¬ 
ter.  For  consistency  with  reported  rate  coefficients,  kinetics  are  detailed  using 
a  temporal  dimension  of  days.  Within  the  CE-QUAL-ICM  code,  kinetics 
sources  and  sinks  are  converted  to  a  dimension  of  seconds  befort  employment 
in  the  mass-conservation  equation. 


Algae 

Sources  and  sinks  of  algae  are: 

Growth  (production) 

Basal  metabolism 

Predation 

Settling 

The  governing  equation  for  algal  biomass  is: 

-  BM  -  PR  -  (6-2) 

H  J 


B  =  algal  biomass,  expressed  as  carbon  (gm  C  m'3) 

P  =  production  (day1) 

BM  =  basal  metabolism  (day1) 

PR  =  predation  (day1) 

WSa  =  algal  settling  velocity  (m  day*1) 

H  =  depth  of  water  column  (m) 

Production 

Production  by  phytoplankton  is  determined  by  the  availability  of  nutrients, 
by  the  intensity  of  light,  and  by  the  ambient  temperature.  The  effects  of  each 
are  considered  to  be  multiplicative: 
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P  =  PMf(N)f(I)f(T) 


(6-3) 


PM  =  production  under  optimal  conditions  (day'1) 
f(N)  =  effect  of  suboptimal  nutrient  concentration  (0  <  f  <  1) 
f(I)  =  effect  of  suboptimal  illumination  (0  <  f  <  1) 
f(T)  =  effect  of  suboptimal  temperature  (0  <  f  <  1) 


Nutrients 

Carbon,  nitrogen,  phosphorus,  and  silica  are  the  primary  nutrients  required 
for  algal  growth.  Inorganic  carbon  is  usually  available  in  excess  and  is  not 
considered.  Silica  is  omitted  since  this  nutrient  is  required  only  by  diatoms 
which  are  not  distinguished  from  other  algae  in  the  model  application.  The 
effects  of  nitrogen  and  phosphorus  on  growth  are  described  by  the  formulation 
commonly  referred  to  as  “Monod  kinetics”  (Monod  1949).  In  the  Monod 
formulation  (Figure  6-1)  growth  is  dependent  upon  nutrient  availability  at  low 
nutrient  concentrations  but  is  independent  of  nutrients  at  high  concentrations. 

A  key  parameter  in  the  formulation  is  the  “half-saturation  concentration”. 
Growth  rate  is  half  the  maximum  when  available  nutrient  concentration  equals 
the  half-saturation  concentration.  Liebig’s  “law  of  the  minimum”  (Odum 
1971)  indicates  growth  is  determined  by  the  nutrient  in  least  supply: 


f(N)  *  minimum 


nh4  +  no3 


PO.d 


KHn  +  NH4  +  N03  KHp  +  P04d 


(6-4) 


NH4  =  ammonium  concentration  (gm  N  m'3) 

N03  =  nitrate  concentration  (gm  N  m'3) 

KHn  =  half-saturation  constant  for  nitrogen  uptake  (gm  N  m  3) 
P04d  =  dissolved  phosphate  concentration  (gm  P  m'3) 

KHp  =  half-saturation  constant  for  phosphorus  uptake  (gm  P  m  3) 


Ught 

Algal  production  increases  as  a  function  of  light  intensity  until  an  optimal 
intensity  is  reached.  Beyond  the  optimal  intensity,  production  declines  as 
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Figure  6-1 .  The  Monod  formulation  for  nutrient-limited  growth 


intensity  increases.  Steele’s  equation  (DiToro  et  al.  1971)  describes  this 
phenomenon: 


f(D 


(b-5) 


I  =  illumination  rate  (Langleys  day'1) 

Is  =  optimal  illumination  (Langleys  day'1) 

Steele’s  equation  describes  the  instantaneous  light  limitation  at  a  point  in 
space.  The  model,  however,  computes  processes  integrated  over  discrete  time 
intervals  and  aggregated  spatially  into  model  segments.  Therefore,  Steele’s 
equation  must  be  integrated  over  an  appropriate  time  interval  and  averaged 
over  the  thickness  of  each  model  segment  The  integration  interval  selected  is 
one  day.  This  interval  does  not  preclude  computation  steps  less  than  a  day  but 
frees  die  model  from  accounting  for  illumination  in  “real  time.”  Daily  averag¬ 
ing  does  preclude  computation  of  diurnal  fluctuations  in  algal  production.  This 
restriction  is  not  severe,  however,  since  the  classic  equations  for  algal  growth 
are  not  appropriate  for  short  time  scales. 
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Assuming  light  intensity  declines  exponentially  with  depth,  the  integrated, 
averaged  form  of  Steele's  equation  is: 


fa> 


2.72  FD  ab 
KessH  ' 


(6-6) 


ab--J2_e-K~H 

FD  Is 


(6-7) 


at  ■  - 


lo 

FDIs 
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Io  =  daily  illumination  at  water  surface  (Langleys  day'1) 

FD  =  fractional  daylength  (0  <  FD  <  1) 

Kess  =  total  light  attenuation  coefficient  (m1) 

Light  attenuation  in  die  water  column  is  composed  ot  two  fractions:  a  back¬ 
ground  value  dependent  on  water  color  and  concentration  of  suspended  parti¬ 
cles,  and  extinction  due  to  light  absorption  by  ambient  chlorophyll: 


Kess  »  Keb  +  Kechl 

CChl 


(6-9) 


Keb  =  background  light  attenuation  (m'1) 

Kechl  =  light  attenuation  coefficient  for  chlorophyll  ’a’  (m2  mg1) 

CChl  =  algal  carbon-to -chlorophyll  ratio  (gm  C  mg'1  chi) 

Optimal  illumination  for  photosynthesis  depends  on  algal  taxonomy,  dura¬ 
tion  of  exposure,  temperature,  nutritional  status,  and  previous  acclimation. 
Variations  in  optimal  illumination  are  largely  due  to  adaptations  by  algae 
intended  to  maximize  production  in  a  variable  environment.  Steele  (1962) 
noted  die  result  of  adaptations  is  that  optimal  illumination  is  a  consistent  frac¬ 
tion  (-  50%)  of  daily  illumination.  Kremer  and  Nixon  (1978)  reported  an 
analogous  finding  that  maximum  algal  production  occurs  at  a  constant  depth 
(-  lm)  in  the  water  column.  Their  approach  is  adopted  here  so  that  optimal 
illumination  is  expressed: 
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-  KessDopt 
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Is  *  Ioavge 


Ioavg  =  adjusted  surface  illumination  (Langleys  day'1) 

Dope  =  depth  of  maximum  algal  production  (m) 

A  minimum,  Ismin,  is  specified  for  optimal  illumination  so  that  algae  do 
not  thrive  at  extremely  low  light  levels.  The  time  required  for  algae  to  adapt 
to  changes  in  illumination  is  recognized  by  computing  Is  based  on  a  time- 
weighted  average  of  daily  illumination: 


Ioavg  =  0.7  Io  ♦  0.2  II  ♦  0.1 12 


(6-11) 


11  =  daily  illumination  one  day  preceding  model  day  (Langleys  day'1) 

12  =  daily  illumination  two  days  preceding  model  day  (Langleys  day'1) 

Four  independent  parameters  influence  the  effect  of  light  on  algal  produc¬ 
tion:  FD,  Ke,  Dope,  and  H.  The  parameters  Ke  and  Dopt,  and  Ke  and  H  occur 
as  products,  however,  so  the  number  of  independent  parameters  that  determine 
the  light  effect  is  actually  three:  FD,  Ke  *  Dopt,  and  Ke  *  H. 

Fractional  daylength,  FD,  occurs  as  a  multiplier  and  in  the  exponential 
terms  of  the  integrated  form  of  Steele’s  equation  (Equation  6-6).  The  net 
effect  of  FD  on  the  light  function  is  nearly  linear,  however  (Figure  6-2).  At 
the  latitude  of  Indian  River,  FD  is  limited  to  the  range  *0.4  to  0.6  so  that  f(I) 
varies  by  «50%  from  the  shortest  to  longest  day  of  the  year. 

A  plot  of  the  effects  on  growth  of  the  pairs  Ke  *  H  and  Ke  ♦  Dopt  (Fig¬ 
ure  6-3)  has  numerous  interpretations.  For  Ke  *  H  »  Ke  *  Dopt,  growth 
diminishes  as  Ke  *  H  increases.  In  other  words,  growth  goes  down  as  light 
extinction  and/or  depth  increase.  For  much  of  the  feasible  range  of  Ke  *  H, 
growth  is  only  weakly  dependent  on  Dopt  This  weak  dependence  was  one 
factor  in  the  selection  of  the  light-effect  formulation.  Maximum  growth  occurs 
when  Ke  *  H  is  slightly  larger  than  Ke  *  Dopt  (when  H  is  slightly  larger  than 
Dopt).  When  Ke  *  H  <  Ke  *  Dopt  (H  <  Dopt),  growth  is  diminished  from  the 
maximum  due  to  light  levels  in  the  majority  of  the  water  column  occupying 
the  supersaturated  range  for  algal  growth.  For  small  Ke  *  H  (shallow  systems 
with  little  light  attenuation)  maximum  growth  occurs  when  Dopt  =  0.  When 
Dopt  =  0,  no  portion  of  the  water  column  is  in  the  supersaturated  range.  For 
much  of  the  feasible  range  of  Ke  *  H,  however,  maximum  production  occurs 
arisen  Dopt  >0.  For  fins  case,  diminished  production  due  to  light  supersatura¬ 
tion  near  the  surface  is  offset  by  increased  production  at  deeper  depths. 
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Figure  6-2.  Effect  of  fractional  dayiength  on  algal  light  limitation 


Figure  6-3.  Effect  of  light  extinction,  total  depth,  and  optimal  depth  on  algal  light  limitation 
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Temperature 

Algal  production  increases  as  a  function  of  temperature  until  an  optimum 
temperature  or  temperature  range  is  reached.  Above  the  optimum,  production 
declines  until  a  temperature  lethal  to  the  organisms  is  attained.  Numerous 
functional  representations  of  temperature  effects  are  available.  Inspection  of 
growth  versus  temperature  curves  indicates  a  function  similar  to  a  Gaussian 
probability  curve  (Figure  6-4)  provides  a  good  fit  to  observations: 


f(T)  =  e  '  KTgl  &  "  Tm)2  when  TSTm 
=  e  ~  KTg2  ^Tm  "  1)2  when  T  >  Tm 


(6-12) 


T  =  temperature  (C°) 

Tm  =  optimal  temperature  for  algal  growth  (C°) 

KTgl  =  effect  of  temperature  below  Tm  on  growth  (C°  "2) 
KTg2  =  effect  of  temperature  above  Tm  on  growth  (C°  2) 


Figure  6-4.  Effect  of  temperature  on  algal  production 
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Basal  Metabolism 


As  employed  here,  basal  metabolism  is  the  sum  of  all  internal  processes 
that  decrease  algal  biomass.  A  portion  of  metabolism  is  respiration  which  may 
be  viewed  as  a  reversal  of  production.  In  respiration,  carbon  and  nutrients  are 
returned  to  the  environment  accompanied  by  the  consumption  of  dissolved 
oxygen.  A  second  internal  sink  of  biomass  is  the  excretion  of  dissolved 
organic  carbon. 

Respiration  cannot  proceed  in  the  absence  of  oxygen.  Basal  metabolism 
cannot  decrease  in  proportion  to  oxygen  availability,  however,  or  algae  would 
approach  immortality  under  anoxic  conditions.  To  solve  this  dilemma,  basal 
metabolism  is  considered  to  be  independent  of  dissolved  oxygen  concentration 
but  the  distribution  of  metabolism  between  respiration  and  excretion  is  oxygen- 
dependent.  When  oxygen  is  freely  available,  respiration  is  a  large  fraction  of 
the  total.  When  oxygen  is  restricted,  excretion  becomes  dominant  Formula¬ 
tion  of  this  process  is  detailed  in  the  text  that  describes  algal  effects  on  carbon 
and  dissolved  oxygen. 

Basal  metabolism  is  commonly  considered  to  be  an  exponentially  increasing 
(Figure  6-5)  function  of  temperature: 


BM  =  BMreKTb(T  "  Tr)  (6-13) 

BMr  =  metabolic  rate  at  Tr  (day*1) 

KTb  =  effect  of  temperature  on  metabolism  (C°  4) 

Tr  =  reference  temperature  for  metabolism  (C°) 


Predation 

Detailed  specification  of  predation  rate  requires  predictive  modeling  of 
zooplankton  biomass  and  activity.  At  present,  zooplankton  are  not  included  in 
the  model.  Consequently,  a  constant  predation  rate  is  specified.  This  specifi¬ 
cation  implicitly  assumes  zooplankton  biomass  is  a  constant  fraction  of  algal 
biomass.  Zooplankton  activity  is  assumed  to  be  influenced  by  temperature  and 
is  taken  into  account  by  incorporating  an  exponential  temperature  relationship 
(Figure  6-5)  into  the  predation  term.  The  predation  formulation  is  identical  to 
basal  metabolism.  The  difference  in  predation  and  basal  metabolism  lies  in  the 
distribution  of  the  end  products  of  these  processes. 


PR  =  BPReK,n,(T  _Tr) 


(6-14) 


BPR  =  predation  rate  at  Tr  (day*1) 
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Figure  6-5.  Exponential  temperature  function 

Effect  of  Algae  on  Organic  Carbon 

During  production  and  respiration,  algae  primarily  take  up  and  produce 
carbon  dioxide,  an  inorganic  form  not  considered  in  the  model.  A  small  frac¬ 
tion  of  basal  metabolism  is  exuded  as  dissolved  organic  carbon,  however,  and 
in  the  model  this  fraction  increases  as  dissolved  oxygen  becomes  scarce. 

Algae  also  produce  organic  carbon  through  the  effects  of  predation.  Zooplank¬ 
ton  take  up  and  redistribute  algal  carbon  through  grazing,  assimilation,  respira¬ 
tion,  and  excretion.  Since  zooplankton  are  not  included  in  the  model,  routing 
of  algal  carbon  through  zooplankton  is  simulated  by  empirical  distribution 
coefficients.  The  effects  of  algae  on  organic  carbon  are  expressed: 


8 

"St 


DOC 


FCD  +  (1  -  FCD) 


KHr 


KHr  +  DO 


\ 

j 


BM  +  FCDPPR 


B 


(6-15) 


J.LPOC  =  FCLPPRB  (6-16) 

^t 
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Ji.  RPOC  =  FCRPPRB 
5t 


(6-17) 


DOC  =  dissolved  organic  carbon  concentration  (gm  C  m  3) 

DO  =  dissolved  oxygen  concentration  (gm  02  m'3) 

LPOC  =  labile  particulate  organic  carbon  concentration  (gm  C  m'3) 

RPOC  =  refractory  particulate  organic  carbon  concentration  (gm  C  m"3) 

FCD  =  fraction  of  basal  metabolism  exuded  as  dissolved  organic  carbon 

KHr  =  half-saturation  concentration  for  algal  dissolved  organic  carbon 
excretion  (gm  02  m‘3) 

FCDP  =  fraction  of  dissolved  organic  carbon  produced  by  predation 
FCLP  =  fraction  of  labile  particulate  carbon  produced  by  predation 
FCRP  =  fraction  of  refractory  particulate  carbon  produced  by  predation 
The  sum  of  the  three  predation  fractions  must  equal  unity. 


Effect  of  Algae  on  Phosphorus 

Algae  take  up  dissolved  phosphate  during  production  and  release  dissolved 
phosphate  and  organic  phosphorus  through  mortality.  As  with  carbon,  the  fate 
of  algal  phosphorus  released  by  metabolism  and  predation  is  represented  by 
distribution  coefficients.  Since  the  total  phosphate  state  variable  includes  both 
intra  and  extracellular  phosphate,  no  explicit  representation  of  the  effect  of 
algae  on  phosphate  is  necessary.  Distribution  of  total  phosphate  is  determined 
by  partition  coefficients  as  detailed  in  the  “PHOSPHORUS”  section  of  this 
chapter.  The  equations  that  express  the  effects  of  algae  on  organic  phosphorus 
are: 


5 

-5? 


DOP  *  (BMFPD  +  PR FPDP)  APC B 


(6-18) 


J.LPOP  *  (BMFPL  +  PR  FPLP)  APC  B 


(6-19) 
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■JLrPOP  =  (BMFPR  ♦  PR FPRP) APC B 


(6-20) 


DOP  =  dissolved  organic  phosphorus  concentration  (gm  P  m'3) 

LPOP  =  labile  particulate  organic  phosphorus  concentration  (gm  P  m  3) 

RPOP  =  refractory  particulate  organic  phosphorus  concentration  (gm  P  m'3) 

APC  =  algal  phosphorus-to-carbon  ratio  (gm  P  gm'1  Q 

FPD  *  fraction  of  dissolved  organic  phosphorus  produced  by  metabolism 

FPL  =  fraction  of  labile  particulate  phosphorus  produced  by  metabolism 

FPR  =  fraction  of  refractory  particulate  phosphorus  produced  by 
metabolism 

FPDP  as  fraction  of  dissolved  organic  phosphorus  produced  by  predation 

FPLP  =  fraction  of  labile  particulate  phosphorus  produced  by  predation 

FPRP  =  fraction  of  refractory  particulate  phosphorus  produced  by  predation 

The  sums  of  the  metabolism  and  respiration  fractions  must  each  be  less 
than  or  equal  to  unity. 


Effect  of  Algae  on  Nitrogen 

Algae  take  up  ammonium  and  nitrate  during  production  and  release  ammo¬ 
nium  and  organic  nitrogen  through  mortality.  Nitrate  is  internally  reduced  to 
ammonium  before  synthesis  into  biomass  occurs  (Parsons  et  aL  1984).  Trace 
concentrations  of  ammonium  inhibit  nitrate  reduction  so  that,  in  the  presence 
of  ammonium  and  nitrate,  ammonium  is  utilized  first  The  “preference”  of 
algae  for  ammonium  can  be  expressed  empirically  (Thom  aim  and  Fitzpatrick 
1982): 


PN  =  NH4. 
+  NH4 


NO, 


(KHn  ♦  NH^)  (KHn  +  N03) 
KHn 

(NH4  ♦  N03)(KHn  ♦  NO3) 


(6-21) 
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PN  =  algal  preference  for  ammonium  uptake  (0  <  PN  <  1) 

The  ammonium  preference  function  (Figure  6-6)  has  two  limiting  values. 
When  nitrate  is  absent,  the  preference  for  ammonium  is  unity.  When  ammo¬ 
nium  is  absent,  the  preference  is  zero.  In  the  presence  of  ammonium  and 
nitrate,  the  preference  depends  on  the  abundance  of  both  forms  relative  to  the 
half-saturation  constant  for  nitrogen  uptake.  When  both  ammonium  and  nitrate 
are  abundant,  the  preference  for  ammonium  approaches  unity.  When  ammo¬ 
nium  is  scarce  but  nitrate  is  abundant,  the  preference  decreases  in  magnitude 
and  a  significant  fraction  of  algal  nitrogen  requirement  comes  from  nitrate. 


Figure  6-6.  Ammonium  preference  function 


The  fate  of  algal  nitrogen  released  by  metabolism  and  predation  is  repre¬ 
sented  by  distribution  coefficients.  The  effects  of  algae  on  the  nitrogen  state 
variables  are  expressed: 


J.NH4  *  (BMFNI  +  PRFNIP  -  PNP)ANCB  (6-22) 

ot 
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^N03  -  (PN  -  1)PANCB 


(6-23) 


£  DON  =  (BM  FND  +  PR  FNDP)  ANC  B  (6-24) 

dt 

-AlPON  =  (BMFNL  ♦  PR  FNLP)  ANC  B  (6-25) 

Ot 

^  RPON  =  (BM  FNR  +  PR  FNRP)  ANC  B  (6-26) 

DON  =  dissolved  organic  nitrogen  concentration  (gm  N  m3) 

LPON  =  labile  particulate  organic  nitrogen  concentration  (gm  N  m'3) 

RPON  =  refractory  particulate  organic  nitrogen  concentration  (gm  N  m'3) 

ANC  =  algal  nitrogen- to-carbon  ratio  (gm  N  gm'1  C) 

FNI  =  fraction  of  inorganic  nitrogen  produced  by  metabolism 

FND  =  fraction  of  dissolved  organic  nitrogen  produced  by  metabolism 

FNL  =  fraction  of  labile  particulate  nitrogen  produced  by  metabolism 

FNR  =  fraction  of  refractory  particulate  nitrogen  produced  by  metabolism 

FNIP  =  fraction  of  inorganic  nitrogen  produced  by  predation 

FNDP  =  fraction  of  dissolved  organic  nitrogen  produced  by  predation 

FNLP  =  fraction  of  labile  particulate  nitrogen  produced  by  predation 

FNRP  =  fraction  of  refractory  particulate  nitrogen  produced  by  predation 

The  sums  of  the  metabolism  fractions  and  the  predation  fractions  must  each 
equal  unity. 
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Algal  Stoichiometry 


Algal  biomass  is  quantified  in  units  of  carbon.  In  order  to  express  the 
effects  of  algae  on  nitrogen  and  phosphorus,  the  ratios  of  nitrogen-to-carbon 
and  phosphorus-to-carbon  in  algal  biomass  must  be  specified.  Global  mean 
values  of  these  ratios  are  well  known  (Redfield  et  al.  1966).  Algal  composi¬ 
tion  varies,  however,  especially  as  a  function  of  nutrient  availability.  As  nitro¬ 
gen  and  phosphorus  become  scarce,  algae  adjust  their  composition  so  that 
smaller  quantities  of  these  vital  nutrients  are  required  to  produce  carbonaceous 
biomass  (Droop  1973;  DiToro  1980;  Parsons  et  al.  1984). 

Observations  from  upper  Chesapeake  Bay  were  examined  to  assess  the 
potential  variability  of  algal  stoichiometry  (Cerco  and  Cole  1993).  The  varia¬ 
tion  of  carbon-to-nitrogen  stoichiometry  in  the  upper  Bay  was  small.  As  a 
consequence  of  these  observations,  the  model  formulation  specified  constant 
algal  nitrogen-to-carbon  ratio,  ANC.  Large  variations  in  carbon-to-phosphorus 
ratio  occurred,  however.  To  account  for  this  effect,  a  variable  algal  phospho- 
rus-to-carbon  ratio,  APC,  was  formulated  in  the  model.  Variable  stoichiometry 
was  not  employed  in  the  Indian  River-Rehoboth  Bay  application.  The  formu¬ 
lation  that  allowed  variable  APC  required  distinctive  treatment  of  the  phos¬ 
phate  state  variable,  however. 


Effect  of  Algae  on  Dissolved  Oxygen 

Algae  produce  oxygen  during  photosynthesis  and  consume  oxygen  through 
respiration.  The  quantity  produced  depends  on  the  form  of  nitrogen  utilized 
for  growth.  More  oxygen  is  produced,  per  unit  of  carbon  fixed,  when  nitrate 
is  the  algal  nitrogen  source  than  when  ammonium  is  the  source.  Equations 
describing  algal  uptake  of  carbon  and  nitrogen  and  production  of  dissolved 
oxygen  (Morel  1983)  are: 


106  CO,  +  16  NH/  +  H,PO;  +  106  H,0  -4 

1  4  2  4  (6-27) 

protoplasm  +  106  02  +  15  H* 


106  C02  -*■  16N03"  +  H2P04"  +  122  H20  +  17  H* 

protoplasm  +  138 


(6-28) 
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When  ammonium  is  the  nitrogen  source,  one  mole  oxygen  is  produced  per 
mole  carbon  dioxide  fixed.  When  nitrate  is  the  nitrogen  source,  1.3  moles 
oxygen  are  produced  per  mole  carbon  dioxide  fixed. 

The  equation  that  describes  the  effect  of  algae  on  dissolved  oxygen  in  the 
model  is: 


( 

(1.3  -  0.3  PN)  P 


DO 


KHr  +  DO 


BM 


AOCRB 


J 


(6-29) 


AOCR  =  dissolved  oxygen-to-carbon  ratio  in  respiration  (2.67  gm  02 
gm1  C) 

The  magnitude  of  AOCR  is  derived  from  a  simple  representation  of  the 
respiration  process: 

CH20  +  02  =  C02  +  H20  (6-30) 


The  quantity  (1.3  -  0.3  PN)  is  die  photosynthesis  ratio  and  expresses  the 
molar  quantity  of  oxygen  produced  per  mole  carbon  fixed.  The  photosynthesis 
ratio  approaches  unity  as  the  algal  preference  for  ammonium  approaches  unity. 


Organic  Carbon 

Organic  carbon  undergoes  innumerable  transformations  in  the  water  col¬ 
umn.  The  model  carbon  cycle  (Figure  6-7)  consists  of  die  following  elements: 

Phytoplankton  production 
Phytoplankton  exudation 
Predation  on  phytoplankton 
Dissolution  of  particulate  carbon 
Heterotrophic  respiration 
Denitrification 
Settling 

Algal  production  is  the  primary  carbon  source  although  carbon  also  enters 
the  system  through  external  loading.  Predation  on  algae  releases  particulate 
and  dissolved  organic  carbon  to  the  water  column.  A  fraction  of  the  particu¬ 
late  organic  carbon  undergoes  first-order  dissolution  to  dissolved  organic  car¬ 
bon.  The  remainder  settles  to  the  sediments.  Dissolved  organic  carbon 
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Figure  6-7.  The  model  cartoon  cycle 


produced  by  phytoplankton  exudation,  by  predation,  and  by  dissolution  is 
respired  or  denitrified  at  a  first-order  rate  to  inorganic  carbon.  No  carbon  is 
recycled  from  the  sediments  to  the  water  column. 
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Dissolution  and  Respiration  Rates 

Dissolution  and  respiration  rates  depend  on  the  availability  of  carbonaceous 
substrate  and  on  heterotrophic  activity.  Heterotrophic  activity  and  biomass 
have  been  correlated  with  algal  activity  and  biomass  across  a  wide  range  of 
natural  systems  (Bird  and  Kalff  1984;  Cole  et  al.  1988).  Consequently,  algal 
biomass  can  be  incorporated  into  dissolution  and  respiration  rate  formulations 
as  a  surrogate  for  heterotrophic  activity.  The  correlation  between  algae  and 
heterotrophs  occurs  because  algae  produce  labile  carbon  that  fuels  heterotro¬ 
phic  activity.  Dissolution  and  respiration  processes  do  not  require  the  presence 
of  algae,  however,  and  may  be  fueled  entirely  by  external  carbon  inputs. 
Representation  of  dissolution  and  respiration  in  the  model  allows  specification 
of  algal-dependent  and  algal-independent  rates: 


Kdoc  =  Kdc  +  KdcalgB 


(6-31) 


Kdoc  =  respiration  rate  of  dissolved  organic  carbon  (day1) 

Kdc  =  minimum  respiration  rate  (day1) 

Kdcalg  =  constant  that  relates  respiration  to  algal  biomass  (m3  gm'1 
C  day1) 


Klpoc  =  Klc  +  KlcalgB 


(6-32) 


Klpoc  =  dissolution  rate  of  labile  particulate  organic  carbon  (day1) 

Klc  =  minimum  dissolution  rate  (day'1) 

Klcalg  =  constant  that  relates  dissolution  to  algal  biomass  (m3  gm' 
C  day1) 


Krpoc  =  Krc  +  Krcalg  B 


(6-33) 


Krpoc  =  dissolution  rate  of  refractory  particulate  organic  carbon  (day1) 

Krc  =  minimum  dissolution  rate  (day'1) 

Krcalg  -  constant  that  relates  dissolution  to  algal  biomass  (m3  gm'1 
C  day1) 

An  exponential  function  (Figure  6-5)  relates  dissolution  and  respiration  to 
temperature. 
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Denitrification 


As  oxygen  is  depleted  from  natural  systems,  oxidation  of  organic  matter  is 
effected  by  the  reduction  of  alternate  oxidants  (referred  to  as  “alternate 
electron  acceptors”).  The  sequence  in  which  alternate  acceptors  are  employed 
is  determined  by  the  thermodynamics  of  oxidation-reduction  reactions.  The 
first  substance  reduced  in  tire  absence  of  oxygen  is  nitrate.  A  representation  of 
the  denitrification  reaction  can  be  obtained  by  balancing  standard  half-cell 
redox  reactions  (Stumm  and  Morgan  1981): 


4N03"  +  4 H *  +  5CH20  — >  2N2  +  7H20  +  5C02  (6*34> 


Equation  6-34  describes  the  stoichiometry  of  the  denitrification  reaction. 
The  kinetics  of  the  reaction,  represented  in  the  model,  are  first-order.  The 
dissolved  organic  carbon  respiration  rate,  Kdoc,  is  modified  so  that  significant 
decay  via  denitrification  occurs  only  when  nitrate  is  freely  available  and  dis¬ 
solved  oxygen  is  depleted  (Figure  6-8).  A  parameter  is  included  so  that  the 
anoxic  respiration  rate  is  slower  than  oxic  respiration: 

Denit  = _ _ _ _  AANOX  Kdoc  (6-35) 

KHodoc  +  DO  KHndn  +  N03 


Denit  =  denitrification  rate  of  dissolved  organic  carbon  (day1) 

AANOX  =  ratio  of  denitrification  to  oxic  carbon  respiration  rate 
(0<  AANOX <  1) 

KHodoc  =  half-saturation  concentration  of  dissolved  oxygen  required  for 
oxic  respiration  (gm  02  tn'3) 

KHndn  =  half-saturation  concentration  of  nitrate  required  for 
denitrification  (gm  N  m'3) 

An  exponential  function  (Figure  6-5)  relates  denitrification  to  temperature. 
Coefficient  values  in  the  function  are  the  same  as  for  dissolved  organic  carbon 
respiration. 

Dissolved  Organic  Carbon 

The  complete  representation  of  all  dissolved  organic  carbon  sources  and 
sinks  in  die  model  ecosystem  is: 
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Denit  /  (AANOX  • 


Figure  6-8.  Effect  of  nitrate  and  dissolved  oxygen  on  denitrification  rate 


Labile  Particulate  Organic  Carbon 

Ttie  complete  representation  of  all  labile  paniculate  organic  carbon  sources 
and  sinks  in  die  model  ecosystem  is: 
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J.LPOC  »  FCLPPRB  -  KlpocLPOC  -  J^LPOC  (6-37) 
Si  H 


WSl  =  settling  velocity  of  labile  particles  (m  day'1) 


Refractory  Particulate  Organic  Carbon 

The  complete  representation  of  all  refractory  particulate  organic  carbon 
sources  and  sinks  in  the  model  ecosystem  is: 


J.RPOC 


FCRPPRB  -  KrpocRPOC  - 


WSr 

H 


RPOC 


(6-38) 


WSr  =  settling  velocity  of  labile  panicles  (m  day'1) 


Phosphorus 

The  model  phosphorus  cycle  (Figure  6-9)  includes  the  following  processes: 

Algal  production  and  metabolism 
Predation 

Hydrolysis  of  particulate  organic  phosphorus 
Mineralization  of  dissolved  organic  phosphorus 
Settling 

External  loads  provide  the  ultimate  source  of  phosphorus  to  die  system. 
Dissolved  phosphate  is  incorporated  by  algae  during  growth  and  released  as 
phosphate  and  organic  phosphorus  through  respiration  and  predation.  A  por¬ 
tion  of  the  particulate  organic  phosphorus  hydrolyzes  to  dissolved  organic 
phosphorus.  Dissolved  organic  phosphorus  is  mineralized  to  phosphate.  Par¬ 
ticulate  organic  phosphorus  which  is  not  mineralized  in  die  water  column 
settles  to  the  sediments.  The  sediments  recycle  phosphorus  to  the  water  col¬ 
umn  as  dissolved  phosphate.  Direct,  mechanistic  linkage  of  phosphorus  release 
to  settling  of  organic  phosphorus  is  not  activated  in  this  model  application, 
however. 

Effects  on  phosphorus  of  algal  production,  metabolism,  and  predation  have 
already  been  detailed.  Descriptions  of  hydrolysis  and  mineralization  and  of  the 
total  phosphate  system  follow. 
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EXTERNAL 

LOADS 


Figure  6-9.  The  mode!  phosphorus  cycle 

Hydrolysis  and  Mineralization 

Within  die  model,  hydrolysis  is  defined  as  the  process  by  which  paniculate 
organic  substances  are  converted  to  dissolved  organic  form.  Mineralization  is 
defined  as  the  process  by  which  dissolved  organic  substances  are  converted  to 
dissolved  inorganic  form.  Conversion  of  particulate  organic  phosphorus  to 
phosphate  proceeds  through  the  sequence  of  hydrolysis  and  mineralization. 
Direct  mineralization  of  particulate  organic  phosphorus  does  not  occur. 

Mineralization  of  organic  phosphorus  is  mediated  by  the  release  of  nucleo¬ 
tidase  and  phosphatase  enzymes  by  bacteria  (Ammerman  and  Azam  1985; 
Chrost  and  Overbeck  1987)  and  algae  (Matavulj  and  Flint  1987;  Chrost  and 
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Overbeck  1987;  Boni  et  al.  1989).  Since  the  algae  themselves  release  the 
enzyme  and  since  bacterial  abundance  is  related  to  algal  biomass,  the  rate  of 
organic  phosphorus  mineralization  is  related,  in  the  model,  to  algal  biomass. 

A  most  remarkable  property  of  the  enzyme  process  is  that  alkaline  phosphatase 
activity  is  inversely  proportional  to  ambient  phosphate  concentration  (Chrost 
and  Overbeck  1987;  Boni  et  al.  1989).  Put  in  different  terms,  when  phosphate 
is  scarce,  algae  stimulate  production  of  an  enzyme  that  mineralizes  organic 
phosphorus  to  phosphate.  This  phenomenon  is  simulated  by  relating  mineral¬ 
ization  to  the  algal  phosphorus  nutrient  limitation.  Mineralization  is  highest 
when  algae  are  strongly  phosphorus  limited  and  is  least  when  no  limitation 
occurs. 

Expressions  for  mineralization  and  hydrolysis  rates  are: 


Kdop  =  Kdp  + _ E!£ _ Kdpalg  B 

KHp  +  P04d 


(6-39) 


Kdop  =  mineralization  rate  of  dissolved  organic  phosphorus  (day1) 

Kdp  =  minimum  mineralization  rate  (day1) 

Kdpalg  =  constant  that  relates  mineralization  to  algal  biomass 
(m3  gm'1  C  day'1) 

Klpop  =  Kip  ♦  KlpalgB  (6-40) 

KHp  +  P04d 


Klpop  =  hydrolysis  rate  of  labile  particulate  phosphorus  (day1) 

Kip  =  minimum  hydrolysis  rate  (day'1) 

Klpalg  =  constant  that  relates  hydrolysis  to  algal  biomass  (m3  gm'1 
C  day1) 

Kipop  =  Kip  Kipalg  B  (Ml) 

KHp  +  P04d 


Kipop  =  hydrolysis  rate  of  refractory  paniculate  phosphorus  (day*1) 

Krp  =  minimum  hydrolysis  rate  (day'1) 

Krpalg  =  constant  that  relates  hydrolysis  to  algal  biomass  (m3  gm'1 
C  day'1) 
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An  exponential  function  (Figure  6-5)  relates  mineralization  and  hydrolysis 
rates  to  temperature. 

Potential  effects  of  algal  biomass  and  nutrient  limitation  on  mineralization 
and  hydrolysis  rates  shown  in  Figure  6-10.  When  nutrient  concentration 
greatly  exceeds  the  half-saturation  concentration  for  algal  uptake,  the  rate 
roughly  equals  the  minimum.  Algal  biomass  has  little  influence.  As  nutrient 
becomes  scarce  relative  to  the  half-saturation  concentration,  the  rate  increases. 
The  magnitude  of  increase  depends  on  algal  biomass.  Factor  of  two  to  three 
increases  are  feasible. 


Figure  6-10.  Effect  of  algal  biomass  and  nutrient  concentration  on  hydrolysis  and 
mineralization 


The  Total  Phosphate  System 

The  model  phosphate  state  variable  is  defined  as  the  sum  of  dissolved  phos¬ 
phate  and  algal  phosphorus  content: 
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P04t  *  P04d  +  P04a 


(6-42) 


P04t  =  total  phosphate  (gm  P  m'3) 

P04d  =  dissolved  phosphate  (gm  P  m'3) 

P04a  =  algal  phosphorus  (gm  P  m  3) 

The  partitioning  of  phosphate  into  extra-cellular  and  intra-cell ular  fractions 
is  required  in  the  formulation  of  variable  phosphorus  stoichiometry.  When 
variable  stoichiometry  is  activated,  algal  phosphate  is  computed  as  a  function 
of  available  nutrient.  In  the  Indian  River-Rehoboth  Bay  application,  however, 
constant  phosphorus  stoichiometry  is  employed.  Algal  phosphate  is  computed: 

P04a  -  APCB  (6-43) 


Dissolved  phosphate  is  the  fraction  of  total  phosphate  not  incorporated  in 
algal  biomass: 


P04d  *  P04t  -  APCB 


(6-44) 


Phosphate 

The  balance  of  the  equations  describing  phosphorus  are  straightforward 
summations  of  previously-described  sources  and  sinks: 

*P04t  -  -  Zp^APCB  +  KdopDOP  (6-45) 

ot  H 


Algal  uptake  and  release  of  phosphate  represents  an  exchange  of  phosphate 
fractions  rather  than  a  phosphate  source  or  sink.  Consequently,  no  algal  source 
or  sink  terms  are  included  in  the  phosphate  mass -conservation  equation.  The 
settling  term  represents  die  settling  of  phosphate  incorporated  in  algal  biomass. 
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Dissolved  Organic  Phosphorus 


^.DOP  =  (BMFPD  +  PR FPDP)  APC B  ♦  KlpopLPOP 

♦  KrpopRPOP  -  KdopDOP 


(6-46) 


Labile  Particulate  Organic  Phosphorus 


5 

Tt 


LPOP 


(BMFPL  ♦  PR FPLP)  APC B  -  KlpopLPOP 


(6-47) 


WS1 


H 


LPOP 


Refractory  Particulate  Organic  Phosphorus 


J.RPOP 

& 


(BMFPR  +  PR FPRP)  APC B  -  KrpopRPOP 


WSr 

H 


RPOP 


(6-48) 


Nitrogen 

The  model  nitrogen  cycle  (Figure  6-1 1)  includes  the  following  processes: 

Algal  production  and  metabolism 
Predation 

Hydrolysis  of  particulate  organic  nitrogen 

Mineralization  of  dissolved  organic  nitrogen 

Settling 

Nitrification 

Denitrification 

External  loads  provide  the  ultimate  source  of  nitrogen  to  the  system.  Inor¬ 
ganic  nitrogen  is  incorporated  by  algae  during  growth  and  released  as 
ammonium  and  organic  nitrogen  through  respiration  and  predation.  A  portion 
of  die  particulate  organic  nitrogen  hydrolyzes  to  dissolved  organic  nitrogen. 


6-28 


Chapter  6  Water  Quality  Modal  Knodca 


Figure  6-1 1 .  The  model  nitrogen  cycle 


The  balance  settles  to  die  sediments.  Dissolved  organic  nitrogen  is  mineralized 
to  ammonium.  In  an  oxygenated  water  column,  a  fraction  of  the  ammonium  is 
subsequently  oxidized  to  nitrate  through  the  nitrification  process.  In  anoxic 
water,  nitrate  is  lost  to  nitrogen  gas  through  denitrification.  Particulate  nitro¬ 
gen  that  settles  to  the  sediments  is  mineralized  and  recycled  to  die  water 
column,  primarily  as  ammonium.  Direct,  mechanistic  linkage  of  ammonium 
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release  to  nitrogen  settling  is  not  activated  in  this  model  application,  however. 
Nitrate  moves  in  both  directions  across  the  sediment-water  interface,  depending 
on  relative  concentrations  in  the  water  column  and  sediment  interstices. 

Effects  on  nitrogen  of  algal  production,  metabolism,  and  predation  have 
already  been  detailed.  Descriptions  of  hydrolysis,  mineralization,  nitrification 
and  denitrification  follow. 


Hydrolysis  and  Mineralization 

In  the  model,  particulate  organic  nitrogen  is  converted  to  the  dissolved 
organic  form  via  hydrolysis.  Dissolved  organic  nitrogen  is  converted  to 
ammonium  through  mineralization.  Conversion  of  particulate  nitrogen  to 
ammonium  proceeds  through  the  sequence  of  hydrolysis  and  mineralization. 
Direct  mineralization  of  particulate  nitrogen  does  not  occur.  The  argument  for 
accelerated  hydrolysis  and  mineralization  during  nutrient-limited  conditions  is 
not  as  clear  for  nitrogen  as  for  phosphorus.  The  same  formulations  are  made 
available  for  nitrogen  as  for  phosphorus,  however.  Accelerated  processes  can 
be  activated  or  deactivated  through  parameter  selection.  The  nitrogen  hydroly¬ 
sis  and  mineralization  formulations  are: 


Kdon  =  Kdn  + _ _ _ Kdnalg  B  (6-49) 

KHn  +  NH4  ♦  N03 


Kdon  =  mineralization  rate  of  dissolved  organic  nitrogen  (day1) 
Kdn  =  minimum  mineralization  rate  (day1) 


Kdnalg  =  constant  that  relates  mineralization  to  algal  biomass 
(m3  gm'1  C  day1) 


Klpon  =  Kin  + 


KHn 

KHn  ♦  NH4  N03 


KlnalgB 


(6-50) 


Klpon  =  hydrolysis  rate  of  labile  particulate  nitrogen  (day1) 

Kin  =  minimum  hydrolysis  rate  (day'1) 

Klnalg  =  constant  that  relates  hydrolysis  to  algal  biomass  (m3  gm'1 

C  day'1) 
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Krpon  =  Km  ♦ 


KHn 


(6-51) 


_ KmalgB 

KHn  ♦  NH4  ♦  N03 


Krpon  =  hydrolysis  rate  of  refractory  particulate  nitrogen  (day1) 

Km  =  minimum  hydrolysis  rate  (day1) 

Kmalg  =  constant  that  relates  hydrolysis  to  algal  biomass  (m3  gm'1 
C  day1) 

An  exponential  function  (Figure  6-5)  relates  mineralization  and  hydrolysis 
rates  to  temperature. 


Nitrification 

Nitrification  is  a  process  mediated  by  specialized  groups  of  autotrophic 
bacteria  that  obtain  energy  through  the  oxidation  of  ammonium  to  nitrite  and 
oxidation  of  nitrite  to  nitrate.  A  simplified  expression  for  complete  nitrifica¬ 
tion  (Tchobanoglous  and  Schroeder  1987)  is: 


NH4  +  2  02  — >  N03"  ♦  H20  «•  2  H  * 


(6-52) 


The  equation  indicates  that  two  moles  of  oxygen  are  required  to  nitrify  one 
mole  of  ammonium  into  nitrate.  The  simplified  equation  is  not  strictly  true, 
however.  Cell  synthesis  by  nitrifying  bacteria  is  accomplished  by  the  fixation 
of  carbon  dioxide  so  that  less  than  two  moles  of  oxygen  are  consumed  per 
mole  ammonium  utilized  (Wezemak  and  Gannon  1968). 

The  kinetics  of  complete  nitrification  are  modelled  as  a  function  of  avail¬ 
able  ammonium,  dissolved  oxygen,  and  temperature: 


NT  = 


DO  NHa 

KHont  “  DO  KHnnt  +  NH4 


f(T)NTm 


(6-53) 


NT  =  nitrification  rate  (gm  N  m'3  day'1) 

KHont  =  half-saturation  constant  of  dissolved  oxygen  required  for 
nitrification  (gm  02  m'3) 
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KHnnt  =  half-saturation  constant  of  NH4  required  for  nitrification 
(gm  N  m'3) 

NTm  =  maximum  nitrification  rate  at  optimal  temperature 
(gm  N  m'3  day*1) 

The  kinetics  formulation  (Figure  6-12)  incorporates  the  products  of  two 
“Monod”  functions.  The  first  function  diminishes  nitrification  at  low  dissolved 
oxygen  concentration.  The  second  function  expresses  the  influence  of  ammo¬ 
nium  concentration  on  nitrification.  When  ammonium  concentration  is  low, 
relative  to  KHnnt,  nitrification  is  proportional  to  ammonium  concentration. 

For  NH4  «  KHnnt,  the  reaction  is  approximately  first-order.  (The  first-order 
decay  constant  =  NTm/KHnnt.)  When  ammonium  concentration  is  large, 
relative  to  KHnnt,  nitrification  approaches  a  maximum  rate.  This  formulation 
is  based  on  a  concept  proposed  by  Tuffey  et  al.  (1974).  Nitrifying  bacteria 
adhere  to  benthic  or  suspended  sediments.  When  ammonium  is  scarce,  vacant 
surfaces  suitable  for  nitrifying  bacteria  exist  As  ammonium  concentration 
increases,  bacterial  biomass  increases,  vacant  surfaces  are  occupied,  and  the 
rate  of  nitrification  increases.  The  bacterial  population  attains  maximum  den¬ 
sity  when  all  surfaces  suitable  for  bacteria  are  occupied.  At  this  point  nitrifi¬ 
cation  proceeds  at  a  maximum  rate  independent  of  additional  increase  in 
ammonium  concentration. 


Figure  6-12.  Effect  of  dissolved  oxygen  and  ammonium  concentration  on  nitrification  rate 
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The  optimal  temperature  for  nitrification  may  be  less  than  peak  tempera¬ 
tures  that  occur  in  coastal  waters.  To  allow  for  a  decrease  in  nitrification  at 
superoptimal  temperature,  the  effect  of  temperature  on  nitrification  is  modelled 
in  the  Gaussian  form  of  Equation  6-4. 


Effect  of  Nitrification  on  Ammonium 


-*LnH4  =  -NT  (6-54) 


Effect  of  Nitrification  on  Nitrate 


-J.NO,  *  NT 
1*  3 


(6-55) 


Effect  of  Nitrification  on  Dissolved  Oxygen 

AdO  *  -  AONTNT  (6-56) 

Ot 


AONT  =  mass  dissolved  oxygen  consumed  per  mass  ammonium-nitrogen 
nitrified  (4.33  gm  O2  gm'1  N) 


Effect  of  Denitrification  on  Nitrate 

The  effect  of  denitrification  on  dissolved  organic  carbon  has  been 
described.  Denitrification  removes  nitrate  from  the  system  in  stoichiometric 
proportion  to  carbon  removal  as  determined  by  Equation  6-34: 

A  NO,  =  -  ANDC  Denit  DOC  (6-57) 

Tit  3 


ANDC  =  mass  nitrate-nitrogen  reduced  per  mass  dissolved  organic  carbon 
oxidized  (0.933  gm  N  gm"1  Q 
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Nitrogen  Mass  Balance  Equations 

The  mass-balance  equations  for  nitrogen  state  variables  are  written  by  sum 
ming  all  previously-described  sources  and  sinks: 

Ammonium 

-*LnH4  =  (BMFNI  +  PRFNIP  -  PNP)ANCB 

(6-58) 

+  Kdon  DON  -  NT 


Dissolved  Organic  Nitrogen 


^LDON  =  (BMFND  +  PR  FNDP)  ANC  B  +  KlponLPON 

+  KrponRPON  -  Kdon  DON 


Labile  Particulate  Organic  Nitrogen 

.AlPON  *  (BMFNL  +  PR FNLP)  ANC B  -  KlponLPON 

(6-60) 

-ZElpon 

H 


Refractory  Particulate  Organic  Nitrogen 


S 

"St 


RPON  =  (BMFNR  +  PRFNRP)  ANC B  -  KrponRPON 


(6-61) 


WSr 

H 


RPON 
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Nitrate 


-*Ln03  =  (PN  -  1)PANCB  +  NT  -  ANDC  Denit  DOC  (6-62) 


Dissolved  Oxygen 

Sources  and  sinks  of  dissolved  oxygen  in  the  water  column  (Figure  6-13) 
include: 

Algal  photosynthesis 
Atmospheric  reaeration 
Algal  respiration 
Heterotrophic  respiration 
Nitrification 


Reaeratlon 

The  rate  of  re  aeration  is  proportional  to  the  dissolved  oxygen  deficit  in 
surface  waters: 


Ado  =  Ai (DOs  -  DO)  (6-63) 

ot  H 

Kr  =  reaeration  coefficient  (m  day'1) 

DOs  =  dissolved  oxygen  saturation  concentration  (gm  02  m'3) 

The  surface  renewal  concept,  attributed  to  Danckwerts  by  O'Connor  and 
Dobbins  (1958),  indicates: 

Kr  =  VD1R  (6_64) 


D1  =  molecular  diffiisivity  of  oxygen  in  water  (*1.7  x  10-4  m2  day'2) 

R  =  surface  renewal  rate  (day'1) 

Specification  of  the  surface  renewal  rate  is  die  fundamental  problem  in 
reaeration  theory.  O'Connor  and  Dobbins  (1958)  state  that,  in  isotropic  turbu¬ 
lence,  surface  renewal  can  be  approximated  as  the  ratio  of  stream  velocity  to 
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ATMOSPHERIC 

OXYGEN 


Figure  6*13.  Dissolved  oxygen  sources  and  sinks 


depth.  The  renewal  rale  is  also  influenced  by  wind,  however  (O’Connor 
1983).  Influences  on  reaeration  of  temperature  (ASCE  1961)  and  salinity 
(Wen  et  al.  1984),  most  likely  effected  through  changes  in  diffusivity,  have 
been  measured.  No  single  theory  that  unites  all  these  factors  into  a  formula¬ 
tion  of  reaeration  in  an  estuary  is  available.  The  surface  renewal  concept  is 
retained  in  this  study  with  the  renewal  rate  treated  as  a  calibration  parameter. 
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* 

Saturation  dissolved  oxygen  concentration  diminishes  as  temperature  and 
salinity  increase.  An  empirical  formula  that  describes  these  effects  (Genet  et 
al.  1974)  is: 


DOs  -  14.5532  -  0.38217  T  ♦  0.0054258  T 2 
-  CL (1.665  xlO"4  -  5.866  xlO^T  ♦  9.796  x  lO^T2) 


(6-65) 


CL  s  chloride  concentration  (=  salinity/1.80655) 

Summary  of  Dissolved  Oxygen  Sources  and  Sinks 

The  complete  kinetics  for  dissolved  oxygen  are: 


"St 


DO 


DO 


(1.3  -  0.3  PN)  P  - _  _  BM 

KHr  +  DO 


-  AONTNT  - 


DO 


KHodoc  +  DO 


AOCRB 

A  OCR  Kdoc  DOC 


(6-66) 


♦  2^1  (DOs  -  DO) 


Salinity 

No  internal  sources  or  sinks  of  salinity  exist 


Temperature 

A  conservation  of  internal  energy  equation  can  be  written  analogous  to  the 
conservation  of  mass  equation.  The  only  source  or  sink  of  internal  energy 
considered  is  exchange  with  the  atmosphere.  Although  solar  radiation  can 
penetrate  several  meters  into  the  water  column,  radiation-induced  increases  in 
internal  energy  are  here  assigned  entirely  to  the  surface  model  layer. 

For  practical  purposes,  die  internal -energy  equation  can  be  written  as  a 
conservation  of  temperature  equation.  Change  of  temperature  due  to  atmo¬ 
spheric  exchange  is  considered  proportional  to  die  temperature  difference 
between  the  water  surface  and  a  theoretical  equilibrium  temperature  (Edinger 
et  al.  1974): 
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(6-67) 


KT 

pCpH 


(Tc  -T) 


Te  =  equilibrium  temperature  (C°) 

KT  =  Heat  exchange  coefficient  (watt  m'2  C° _1) 

Cp  =  specific  heat  of  water  (4200  watt  sec  kg'1  C°  *1) 
p  =  density  of  water  (1000  kg  m'3) 


Summary  of  Kinetics  Coefficients 

Initial  values  of  kinetics  coefficients  were  obtained  from  the  Chesapeake 
Bay  model  (Cerco  and  Cole,  1993).  These  had  been  determined  following  an 
exhaustive  literature  search  and  extensive  calibration  of  the  model.  Many  of 
the  initial  values  were  suited  for  Indian  River-Rehoboth  Bay  with  no  modifica¬ 
tion.  Values  of  remaining  coefficients  were  determined  through  a  recursive 
process  in  which  initial  values  were  altered  until  optimal  agreement  between 
observations  and  model  predictions  was  achieved.  Coefficient  values  are  sum¬ 
marized  in  Table  6-2.  The  table  also  indicates  whether  coefficients  were 
adopted  from  Chesapeake  Bay  (CB)  or  determined  specifically  for  Indian 
River-Rehoboth  Bay  (IRRB). 
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Table  6-2 

Kinetics  Coefficients 

n 

Symbol 

Definition 

Value 

Origin 

AANOX 

ratio  of  denitrification  to  oxic  carbon 
respiration  rate 

0.5 

CB 

ANC 

algal  nitrogen-to-cartoon  ratio 

0.167  gm  N  gm'1  C 

CB 

APC 

algal  phosphorus-to-carbon  ratio 

0.0167  gm  P  gm'1  C 

IRRB 

BMr 

basal  metabolic  rate  at  Tr 

0.1  day'1 

IRRB 

BPR 

predation  rale  at  Tr 

0  .125  day  '1 

IRRB 

CChl 

algal  carbon-to-chlorophyB  ratio 

60  gm  C  mg'1  chi 

CB 

Oopl 

depth  of  maximum  algal  production 

0.5  m 

IRRB 

FCD 

fraction  of  basal  metabolism  exuded  as 
dissolved  organic  carbon 

0.0 

CB 

FCDP 

fraction  of  dissolved  organic  carbon 
produced  by  predation 

0.0 

IRRB 

FCLP 

fraction  of  labile  particulate  carbon 
produced  by  predation 

0.7 

IRRB 

FCRP 

fraction  of  refractory  particulate  carbon 
produced  by  predation 

0.3 

IRRB 

FNI 

fraction  of  inorganic  nitrogen  produced 
by  metabofam 

0.2 

IRRB 

FNIP 

fraction  of  inorganic  nitrogen  produced 
by  predation 

0.0 

IRRB 

FND 

fraction  of  dissolved  organic  nitrogen 
produced  by  metabolism 

0.8 

IRRB 

FNDP 

fraction  of  dissolved  organic  nitrogen 
produced  by  predation 

0.0 

IRRB 

FNL 

fraction  of  labile  particulate  nitrogen 
produced  by  metaboism 

0.0 

IRRB 

FNLP 

fraction  of  labile  particulate  nitrogen 
produced  by  predation 

0.7 

IRRB 

FNR 

fraction  of  refractory  particulate  nitrogen 
produced  by  metabolism 

0.0 

IRRB 

FNRP 

fraction  of  refractory  particulate  nitrogen 
produced  by  predation 

0.3 

IRRB 

FPD 

fraction  of  dissolved  organic  phosphorus 
produced  by  metabolism 

0.8 

IRRB 

FPDP 

fraction  of  dissolved  organic  phosphorus 
produced  by  predation 

0.0 

IRRB 

FPL 

fraction  of  labile  partictiate  phosphorus 
produced  by  metabotism 

0.0 

IRRB 
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1  Table  6*2  (Continued) 

Symbol 

Definition 

Value 

Origin 

FPU* 

fraction  of  labile  particUale  phosphorus 
produosd  by  predation 

0.7 

IRRB 

FPR 

fraction  at  refractory  particulate 
phosphorus  produced  by  matabotam 

0.0 

IRRB 

FPRP 

fraction  at  refractory  particUate 
phosphorus  produced  by  predation 

0.3 

IRRB 

tamn 

minimum  ilumination  tor  algal  growth 

40  Langleys  day*1 

CB 

Kdc 

minimum  respiration  rale  of  dissolved 
organic  carbon 

0.1  day1 

IRRB 

Kdcalg 

constant  that  rotates  respiration  rale  to 
algal  biomass 

0.0  m3  gm'1  C  ttey'1 

CB 

Kdn 

minimum  mineraitzalion  rate  of  dssofced 
organic  nitrogen 

0.1  day1 

IRRB 

Kdnrig 

constant  that  relates  mineraMzation  rate 
to  algal  biomass 

0.0  m3  gm"1  C  day'1 

CB 

Kdp 

minimum  mineralization  rate  of  dasofved 
organic  phosphorus 

0.1  day'1 

CB 

Kdp eg 

constant  that  relates  mineralization  rate 
to  algal  bioman 

0.2  m3  gm'1  C  thy'1 

CB 

Kechl 

light  attenuation  coefficient  for 
chlorophyll  'a' 

16  rr?  mg'1 

IRRB 

KHn 

half-saturation  concentration  for  nitrogen 
uptake 

0.01  gm  N  m*3 

CB 

KHndn 

half-aaturation  concentration  of  nitrate 
required  tor  denitrification 

0.1  gm  N  m*3 

CB 

KHrmt 

half-saturation  concentration  of  NH4 
required  tor  nitrificalion 

1.0  gm  N  m"3 

CB 

KHodoc 

half-saturation  concentration  of  dissolved 
oxygen  required  for  oxic  rsspiratnn 

0.5  gmO^m’3 

CB 

KHont 

half-saturation  concentetion  of  dissolved 
oxygen  raQuivsd  for  nitrification 

I.OgmO^m"3 

CB 

KHp 

half-saturation  eonoentration  for 
phosphorus  uptake 

0.001  gm  P  m-3 

CB 

KHr 

half -saturation  concentration  for  aigal 
dissolved  organic  carbon  excretion 

0.5  gm  m*3 

CB 

1 

minimum  dbsofution  rate  of  labile 

0.075  day'1 

CB 

I  particulate  carbon 
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(  Table  6-2  (Continued) 

1  Symbol 

Definition 

Value 

Origin 

Kicaig 

constant  that  relates  dissolution  rate  to 
algal  biomass 

0.0  m3  gm'1  C  day'1 

CB 

Kin 

minimum  dissolution  rate  of  labile 
particulate  nitrogen 

0.075  day'1 

CB 

Klnalg 

constant  that  relates  dissolution  rate  to 
algal  biomass 

0.0  m3  gm'1  C  day'1 

CB 

Kip 

minimum  dissolution  rate  of  labile 
particulate  phosphorus 

0.075  day'1 

CB 

KlpaJg 

constant  that  relates  dissolution  rate  to 
algal  biomass 

0.0  m3  gm'1  C  day'1 

CB 

Kr 

reaeration  coefficient 

0.63  m  day'1 

IRRB 

Krc 

minimum  dissolution  rate  of  refractory 
particulate  carbon 

0.005  day'1 

CB 

Krcalg 

constant  that  relates  dissolution  rate  to 
algal  biomass 

0.0  m3  gm'1  C  day'1 

CB 

Km 

minimum  dissolution  rate  of  refractory 
particulate  nitrogen 

0.005  day'1 

CB 

KmaJg 

constant  that  relates  dissolution  rate  to 
algal  biomass 

0.0  m3  gm'1  C  day'1 

CB 

Krp 

minimum  dissolution  rate  of  refractory 
particulate  phosphorus 

0.005  day1 

CB 

Krpalg 

constant  that  relates  dissolution  rate  to 
algal  biomass 

0.0  m3  gm'1  C  day'1 

CB 

KTb 

effect  of  temperature  on  basal 
metabolism 

0.068  C°  *’ 

CB 

KTfll 

effect  of  temperature  below  Tm  on 
growth 

0.012  C0"2 

IRRB 

KTfl2 

effect  of  temperature  above  Tm  on 
growth 

0.012  C0"2 

IRRB 

KThdr 

constant  that  relates  hydrolysis  rates  to 
temperature 

0.068  C°  -1 

CB 

KTmnl 

constant  that  relates  mineralization  rates 
to  temperature 

0.068  C°  0 

CB 

KTntl 

effect  of  temperature  below  Tmnt  on 
nitrification 

0.08  C°  •* 

IRRB 

(Sheets of 4)  | 
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1  Table  6-2  (Concluded) 

||  Symbol 

Definition 

Value 

Origin 

KTnt2 

effect  of  temperature  above  Tmnt  on 
nitrification 

0.09  C0'* 

IRRB 

NTm 

maximum  nitrification  rate  at  optimal 
tamparatura 

0.1  gm  N  m"3  day'1 

IRRB 

PM 

production  under  optimal  conditions 

2.5  day'1 

CB 

Tm 

optimal  tamparatura  for  algal  growth 

25  C° 

CB 

Tmnt 

optimal  tamparatura  for  nitrification 

30  C° 

IRRB 

Tr 

reference  temperature  for  metabolism 

20  C° 

CB 

Trhdr 

reference  temperature  for  hydrolysis 

20  C° 

CB 

Trmnl 

reference  temperature  tor  mineralization 

20  C° 

CB 

WSa 

algal  sealing  velocity 

0.1  mday'1 

CB 

WSI 

sealing  velocity  of  labile  parades 

0.1  m  day'1 

CB 

WSr 

sealing  velocity  of  refractory  particles 

0.1  mday'1 

CB 

I  (Shaal  4  of  4)  | 
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Chapter  VII:  Linking  the 
Hydrodynamic  and  Water 
Quality  Models 


Introduction 

The  hydrodynamic  and  water  quality  models  are  indirectly  linked.  Cell 
dimensions,  flows,  and  other  information  from  the  hydrodynamic  model  are 
written  to  disk  for  subsequent  use  by  the  water  quality  model  The  indirect 
linkage  makes  efficient  use  of  computer  resources.  An  unlimited  number  of 
water  quality  model  runs  can  be  made  using  information  put  out  by  a  single 
run  of  the  hydrodynamic  model.  In  order  to  link  the  models,  output  from  the 
hydrodynamic  model  must  be  ‘‘mapped"  to  die  water  quality  model.  That  is. 
information  in  the  x-y  plane  employed  by  CH3D  must  be  converted  into 
dimension-free  information  for  use  in  die  control  volumes  of  CE-QUAL-ICM. 

The  simplest  mapping  is  one-to-one  in  space  and  time.  In  one-to-one  map¬ 
ping,  water  quality  model  control  volumes  correspond  in  number  and  size  to 
hydrodynamic  grid  cells.  Flows  and  other  information  are  stored  after  each 
time  step  of  die  hydrodynamic  model  and  die  two  models  employ  die  same 
time  step.  One-to-one  mapping  is  the  least  efficient  linkage  procedure.  Vast 
amounts  of  disk  space  are  required  to  store  all  hydrodynamic  information.  The 
water  quality  model  devotes  an  unnecessary  proportion  of  execution  time  to 
reading  hydrodynamic  information  and  integrating  the  conservation  of  mass 
equation  using  the  hydrodynamic  time  step.  For  lengthy  executions  of  the 
hydrodynamic  and  water  quality  models,  the  storage  volume  and  execution 
time  demanded  by  one-to-one  mapping  are  prohibitive. 

Averaging  hydrodynamic  output  before  disk  storage  is  one  technique 
employed  to  increase  linkage  efficiency.  A  second  technique  is  the  grid  over¬ 
lay  in  which  multiple  hydrodynamic  model  cells  are  combined  into  single 
water  quality  model  control  volumes.  Both  averaging  and  overlays  were 
employed  in  the  linkage  applied  to  Indian  River-Rehoboth  Bay.  Hydrodyna¬ 
mic  information,  computed  at  thirty-second  intervals,  was  averaged  into  hourly 
intervals  before  disk  storage.  The  2063  CH3D  cells  were  overlaid  into  a 
CE-QUAL-ICM  grid  that  contained  281  control  volumes  (Figure  7-1). 
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Figure  7-1 .  Water  Quality  Model  Grid 


Water  Quality  Model  Time  Step 

The  interval  at  which  hydrodynamic  output  is  written  is  not  necessarily  the 
appropriate  time  step  for  use  in  the  water  quality  model.  Rather,  the  water 
quality  model  time  step  is  determined  by  the  stability  requirements  of  the 
numerical  scheme  employed  to  integrate  the  conservation  of  mass  equation 
(Equation  6-1).  Typically,  the  time  step  determined  by  stability  limits  is 
shorter  than  the  output  interval.  An  “autostepping”  algorithm  within  the  water 
quality  model  determines  the  time  step.  The  algorithm  computes  permissible 
time  step  based  on  flow,  dispersion,  and  cell  dimension.  If  the  time  step  is 
shorter  than  the  hydrodynamic  output  interval,  flow  and  dispersion  are  held 
constant  until  hydrodynamic  information  is  updated  at  the  next  read  operation. 
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Inspection  of  the  stability  region  of  the  one-dimensional  QUICKEST  algo¬ 
rithm  (Leonard  1979)  indicates  sufficient  conditions  are: 


At  < 


Ax 

T 


and 


At  <  ill 
-  2D 


(7-1) 


(7-2) 


At  =  time  step  (T) 

Ax  =  cell  length  (L) 
u  =  velocity  (L  T1) 

D  =  diffusion  coefficient  (L2  T1) 


The  autostepping  algorithm  examines  velocity,  diffusion,  and  cell  length  at 
each  flow  face  of  die  water  quality  model  control  volumes.  (These  correspond 
to  Qj  /  Aj,  Dj,  and  5xj  in  Equation  6-1).  The  time  step  is  determined  as: 


At  *  minimum 


v  u  2D, 


(7-3) 


a  =  constant  that  insures  time  step  is  less  than  maximum  allowed  («0.9S) 

The  flow  face  with  the  most  restrictive  time  step  determines  the  time  step 
for  the  entire  system. 

Although  the  criteria  expressed  in  Equations  7-1  and  7-2  are  not  strictly 
correct  for  application  of  QUICKEST  in  multiple  dimensions,  we  have  thus  far 
found  the  one-dimensional  criteria  sufficient  to  determine  the  time  step  in 
multi-dimensional  applications.  The  mean  time  step  determined  for  die  water 
quality  model  was  =300  seconds,  an  order  of  magnitude  greater  than  the 
hydrodynamic  model  time  step. 
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Linkage  Testing 


Linkage  between  the  two  models  is  accomplished  through  computer  code 
which  is  not  detailed  here.  Once  the  linkage  software  is  complete,  the  linkage 
must  be  tested  to  insure  that  transport  in  the  water  quality  model  compares  to 
transport  in  the  hydrodynamic  model.  For  one-to-one  linkage  the  testing  is 
primarily  a  test  of  computer  code.  When  hydrodynamic  information  is 
correctly  mapped  into  the  water  quality  model,  transport  in  the  two  models 
should  be  identical.  When  temporal  averaging  and  overlays  are  employed, 
factors  in  addition  to  coding  require  consideration.  As  the  averaging  interval 
and  overlay  size  increase,  transport  in  the  water  quality  model  di  erges  from 
the  hydrodynamic  model.  Judgement  is  required  to  determine  the  acceptable 
temporal  interval  and  spatial  resolution. 

The  ultimate  linkage  test  is  comparison  of  conservative  substance  transport 
in  the  hydrodynamic  and  water  quality  models.  To  perform  our  test,  we  set 
initial  salinity  in  the  hydrodynamic  model  to  zero.  The  salinity  at  the  interface 
with  the  ocean  was  fixed  at  50  ppt  The  hydrodynamic  model  was  run  for 
t’urty  days  during  which  salt  was  advected  and  dispersed  from  the  ocean  inter¬ 
face  the  interior  of  Indian  River  and  Rehoboth  Bay.  We  ran  the  water 
ciual  i ,  vudel  using  flow  and  dispersion  from  the  hydrodynamic  model  and 
identic*  initial  and  boundary  conditions  for  salinity.  One-hour  averaging  of 
hydrodynamic  information  and  the  overlay  grid  were  employed.  Time  series 
of  salinity  predicted  by  the  two  models  were  compared  at  several  locations  on 
the  water  quality  model  grid  (Figure  7-2).  Each  comparison  (Figures  7-3  to 
7-7)  comprised  four  plots.  Time  series  of  salinity  predictions,  output  each 
hour,  were  plotted  in  the  upper  left.  Hourly  salinity  predictions  from  the  two 
models  were  plotted  against  each  other  in  the  upper  right.  Salinity  from  both 
models  was  averaged  over  two  tidal  cycles  (=25  hours)  and  plotted  as  time 
series  in  the  lower  left.  Tidal-average  salinity  predictions  were  plotted  against 
each  other  in  the  lower  right 

Perfect  linkage  is  indicated  when  the  salinity  time  series  from  the  two  mod¬ 
els  exactly  superimpose.  On  the  corresponding  scatterplots,  all  points  will  fall 
on  the  one-to-one  diagonal  line.  The  ideal  is  nearly  realized  at  cells  52,  127, 
and  134.  These  cells  lie  in  regions  in  which  the  horizontal  gradient  of  salinity 
is  small.  At  cells  247  and  1 17,  which  lie  in  regions  of  substantial  horizontal 
salinity  gradient,  agreement  is  less  than  ideal  on  the  hourly  time  scale.  When 
results  are  averaged  over  two  tidal  cycles,  however,  agreement  is  much 
improved. 

The  linkage  tests  indicate  the  temporal  averaging  of  hydrodynamic  output 
and  the  grid  overlay  are  satisfactory  for  transport  problems  in  which  time 
scales  equal  or  exceed  two  tidal  cycles  and  spatial  concentration  gradients  are 
mild.  These  conditions  are  typical  in  eutrophication  studies.  We  conclude  the 
linkage  between  the  two  models  is  correct  and  suited  for  eutrophication  model¬ 
ing  in  Indian  River-Rehoboth  Bay. 
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Figure  7-2.  Cells  Selected  tor  Transport  Comparisons 
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Figure  7-3.  Salinity  Comparison  at  Cell  134 
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Figure  7-4.  SaNnity  Comparison  at  Ceil  127 


ire  7-5.  SaNnfty  Comparison  at  Cell  52 
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Instantaneous  “  W®1  Celt  117 
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Figure  7-7.  Salnlly  Comparison  at  Cell  117 


Chapter  VIII:  Additional  Model 
Inputs  and  Outputs 


Introduction 

Additional  information,  beyond  loads  and  calibration  parameter  values,  is 
required  to  operate  the  model.  Additional  inputs  considered  here  are  load 
partition,  open-mouth  boundary  conditions,  light  extinction,  and  sediment-water 
fluxes.  Postprocessing  of  model  state  variables  into  more  useful  quantities  is 
also  desirable.  Derivation  of  BOD  concentration  and  calculation  of  diurnal 
fluxes  from  model  computations  are  detailed. 


Partition  of  Loads 


The  quantities  analyzed  in  the  free  flowing  streams  and  in  the  point-source 
discharges  required  mapping  into  model  state  variables  (Table  8-1). 


Table  8-1  | 

Correspondence  of  Loads  and  Model  State  Variables  | 

Observation 

maps  into 

Sms  Variable 

Biochemical  Oxygen 
Demand 

- > 

Dissolved  Organic  Carbon 

Labile  Particulate  Organic  Carbon 

Refractory  Particulate  Organic  Carbon 

Total  Kjeldahl  Nitrogen 

- > 

Dissolved  Organic  Nitrogen 

Labile  Particulate  Organic  Nitrogen 

Refractory  Particulate  Organic  Nitrogen 

Total  Phosphorus 

- > 

Dissolved  Organic  Phosphorus 

Labile  Particulate  Organic  Phosphorus 

_ ! _ 

Refractory  Particulate  Organic  Phosphorus 
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No  exact  conversion  of  BOD  to  organic  carbon  was  possible  but  reasonable 
empirical  correspondences  were  employed  (Chapter  3).  The  derived  total 
organic  carbon  still  required  partitioning  into  model  state  variables,  however. 
Relationships  used  to  derive  loads  of  model  state  variables  from  observed 
loads  were: 


DOC  «  FDOC  *  TOC 

(8-1) 

LPOC  *  FLPOC  *  (TOC  -  DOC) 

(8-2) 

RPOC  =  (l  -  FLPOC)  *  (TOC  -  DOC) 

(8-3) 

TOC  =  total  organic  carbon  (gm  m*3) 

DOC  =  dissolved  organic  carbon  (gm  m*3) 

LPOC  =  labile  particulate  organic  carbon  (gm  m'3) 

RPOC  =  refractory  particulate  organic  carbon  (gm  m'3) 

FDOC  =  fraction  of  total  organic  carbon  that  is  dissolved 

FLPOC  =  fraction  of  particulate  organic  carbon  that  is  labile 

DON  *  FDON  *  (TKN  -  NH4) 

(8-4) 

LPON  =  FLPON  ♦  (TKN  -  NH4  -  DON) 

(8-5) 

RPON  =  (l  -  FLPON)  *  (TKN  -  NH4  -  DON) 

(8-6) 

TKN  =  total  Kjeklahl  nitrogen  (gm  m'3) 

NH4  =  ammonium  (gm  m'3) 

DON  =  dissolved  organic  nitrogen  (gm  m'3) 

LPON  =  labile  particulate  organic  nitrogen  (gm  m'3) 
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RPON  *  refractory  particulate  organic  nitrogen  (gm  m'3) 

FDON  =  fraction  of  total  organic  nitrogen  that  is  dissolved 
FLPON  =  fraction  of  particulate  organic  nitrogen  that  is  labile 

DOP  -  FDOP  *  (TOTP  -  P04)  (8'7) 

LPOP  -  FLPOP  *  (tOTP  -  P04  -  DOP)  (8‘8) 

RPOP  -  (l  -  FLPOP)  *  (TOTP  -  P04  -  DOP)  (8'9) 


TOTP  =  total  phosphorus  (gm  m'3) 

P04  =  dissolved  phosphate  (gm  m3) 

DOP  =  dissolved  organic  phosphorus  (gm  m'3) 

LPOP  *  labile  particulate  organic  phosphorus  (gm  m3) 

RPOP  =  refractory  particulate  organic  phosphorus  (gm  m'3) 

FDOP  =  fraction  of  total  organic  phosphorus  that  is  dissolved 

FLPOP  =  fraction  of  particulate  organic  phosphorus  that  is  labile 

Initial  values  of  splits  were  adopted  from  values  employed  in  Chesapeake 
Bay  tributaries.  These  values  were  refined  during  calibration  of  die  Indian 
River-Rehoboth  Bay  model  (Table  8-2). 


Open-Mouth  Boundary  Conditions 

Concentrations  at  the  junction  of  Indian  River  Inlet  with  the  Atlantic  Ocean 
were  specified  primarily  using  data  collected  offshore  of  die  inlet  by  die  EPA 
(Chapter  2).  Observations  were  available  for  12  of  the  36  months  covered  by 
model  calibration.  Boundary  conditions  in  months  for  which  no  observations 
were  available  were  filled  in  with  values  observed  in  die  same  month  in  alter¬ 
nate  years  or  values  observed  at  die  mouth  of  Chesapeake  Bay.  Good  corre¬ 
spondence  was  observed  in  temperature,  salinity,  chlorophyll,  phosphate,  and 
dissolved  oxygen  between  Observations  collected  offshore  of  Indian  River  and 
at  die  bottom  of  the  Bay  mouth  (Figures  8-1  to  8-5).  Salinity  boundary 
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Table  8-2 

Dissolved  and  Labile  Fractions  of  Loads 

Parameter 

Nonpoint-Souice  Value 

Point-Source  Value 

FDOC 

0.8 

0.8 

FLPOC 

0.0 

0.15 

FDON 

0.6 

0.8 

FLPON 

0.0 

0.1 

FDOP 

0.1 

0.4 

FLPOP 

0.0 

0.0 

conditions  for  July  and  August  1989  were  specified  to  match  observations 
collected  within  Indian  River  Bay  during  those  unusually  wet  months. 


Light  Extinction 

Light  extinction  was  specified  based  on  disk  visibility  measures  collected 
by  the  Academy  of  Natural  Sciences  and  the  College  of  Marine  Studies  (Chap¬ 
ter  2).  Disk  visibility  was  converted  to  light  extinction  via  the  relationship: 


(8-10) 


KE  =  light  extinction  (m'1) 

DV  =  secchi  depth  (m) 

Light  extinction  in  Indian  River-Rehoboth  Bay  is  strongly  affected  by  chlo¬ 
rophyll  concentration  (Figure  8-6).  Light  extinction  input  to  the  model  is 
background  extinction  independent  of  chlorophyll.  The  effect  of  chlorophyll  is 
added  to  the  background  value  based  on  computed  chlorophyll  concentration 
(Equation  6-9).  Consequently,  the  effect  of  ambient  chlorophyll  had  to  be 
removed  horn  the  observed  light  extinction.  Light  extinction  in  die  study 
system  may  also  vary  seasonally  or  regionally.  In  order  to  quantify  the  influ¬ 
ence  of  chlorophyll  and  detect  seasonal  and  regional  trends,  a  model  was 
proposed: 
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1  2  3  4  5  6  7  *  9  10  11  12 

MONTH 


Figure  8*1 .  Temperature  at  Mouth  of  Indian  River  and  Chesapeake  Bay 


Figure  8-2.  Salinity  at  Mouth  of  Indian  River  and  Chesapeake  Bay 
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Figure  8-3.  Chlorophyll  at  Mouth  of  Indian  River  and  Chesapeake  Bay 


Figure  8-4.  Phosphate  at  Mouth  of  Indian  River  and  Chesapeake  Bay 
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D0(gm/m**3) 


Figure  8-5.  Dissolved  Oxygen  at  Mouth  of  Indian  River  and  Chesapeake  Bay 


Figure  8-6.  Relation  of  Light  Extinction  to  Chlorophyll  Concentration 
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KE  =  a  +  bRl  ♦  cR2  ♦  dR3 


+  esin 


2  n 


36S 


Jday 


f 


2k 


♦  fco^Jday 


(8-11) 


gCHL 


a  =  annual  mean  background  light  extinction  (m'1) 
b  =  incremental  light  extinction  in  Region  1  (m1) 
c  =  incremental  light  extinction  in  Region  2  (m1) 
d  =  incremental  light  extinction  in  Region  3  (m*1) 

R1  =  1  in  Region  1,  0  elsewhere 
R2  =  1  in  Region  2,  0  elsewhere 
R3  =  1  in  Region  3, 0  elsewhere 

e,  f  =  parameters  that  determine  amplitude  and  phase  of  seasonal  variation 
(m-1) 

Jday  =  julian  day 

g  =  effect  of  chlorophyll  on  light  extinction  (m2  mg’1  Chl’a’) 

Stepwise  regression  was  employed  to  evaluate  parameter  values.  Observed 
extinction  was  treated  as  die  dependent  variate  while  region,  julian  day,  and 
chlorophyll  were  treated  as  independent  variates.  Three  regions  were  defined 
(Figure  8-7):  Region  1  -  Lower  Indian  River,  Region  2  -  Upper  Indian  River, 
Region  3  -  Rehoboth  Bay. 

The  analysis  (Table  8-3)  indicated  that  chlorophyll  was  die  primary  deter¬ 
minant  of  light  extinction  in  Indian  River-Rehoboth  Bay.  A  seasonal  trend, 
independent  of  chlorophyll  was  also  detected.  Background  extinction  in  Upper 
Indian  River  was  greater  than  elsewhere  but  that  no  difference  existed  between 
Lower  Indian  River  and  Rehoboth  Bay.  Monthly  values  of  background  light 
extinction,  that  accounted  for  seasonal  and  regional  variation  (Figure  8-8), 
were  input  to  the  model. 


Sediment-Water  Fluxes 

Sediment-water  fluxes  of  oxygen  and  nutrients  can  be  a  significant  influ¬ 
ence  on  water  quality  in  shallow  systems  such  as  Indian  River-Rehoboth  Bay. 
One  component  of  the  sediment-water  system  in  the  shallow  bays  is  a  layer  of 
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Figure  8-7.  Regions  in  Light  Extinction  Analysis 


Table  8-3 

Analysis  of  Light  Extinction 

Parameter 

Value 

Cumulative  R2 

a 

1.40 

g 

0.016 

0.30 

e 

-056 

0.35 

f 

-0.74 

0.40 

c 

0.88 

0.47 

benthic  algae  that  occupies  the  sediment-water  interface  in  regions  where  light 
intercepts  the  bottom.  Actions  of  the  benthic  algae  modify  the  sediment-water 
fluxes  that  usually  occur  in  a  darkened  water  column.  Photosynthesis  by  ben¬ 
thic  algae  may  induce  an  apparent  sediment  release  of  dissolved  oxygen 
instead  sediment  oxygen  demand.  Nutrient  uptake  by  the  algae  may  induce 
apparent  sediment  nutrient  uptake  when  release  to  the  water  column  is 
expected. 
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Figure  8-8.  Seasonal  and  Regional  Variation  in  Light  Extinction 


A  series  of  empirical  models  of  sediment-water  fluxes  were  constructed 
based  on  observations  collected  in  Indian  River-Rehoboth  Bay  (Chapter  2). 
Raw  data  supplied  by  the  principal  investigator  consisted  of  replicate  measures 
collected  at  varying  light  intensities.  The  replicate  measures  were  averaged 
and  scientific  units  were  converted  to  units  employed  in  the  model.  Simulta¬ 
neous  measures  collected  in  light  and  dark  were  combined  into  daily-average 
fluxes  that  occur  during  day-night  illumination  cycles.  The  observations, 
treated  for  analysis,  are  reported  in  Table  8-4. 

Empirical  models  were  formulated  based  on  examination  of  the  data,  on 
expected  actions  of  benthic  algae,  and  on  typical  sediment-water  interactions  in 
a  darkened  water  column.  The  sediment  oxygen  demand  model  was: 


SOD  =  SODbeoCr  "  20)  ♦  bLight 


(8-12) 


SOD  =  sediment  oxygen  demand  (gm  m*2  day*1) 


SODb  =  sediment  oxygen  demand  at  20  C°  in  a  darkened  water  column 
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Table  8-4 

Observed  Sediment  Water  Fluxes 

Station 

Tamp, 

C° 

Light, 

langtey 

day'1 

SOD, 
gm  m'2 
day'1 

mg  m 
day'1 

NO, 
mg  m*2 
day'1 

mg  m* 
day1 

1 

12 

157.5 

0 

0 

12 

0 

-0.784 

52 

0.45 

2 

12 

157.5 

0.505 

-5.7 

-0.76 

12 

0 

-1.25 

34 

0 

3 

12 

157.5 

0.233 

-5.5 

-0.22 

4 

12 

6.75 

-0.894 

395 

-16 

1 

25 

63 

-1.12 

27 

1.44 

2 

25 

15.3 

-1.06 

34 

856 

25 

0 

-1.54 

ISO 

10.4 

3 

25 

90 

-0.209 

5.9 

0 

25 

0 

-1.15 

107 

6.45 

4 

25 

0 

-3.48 

147 

9.67 

o  =  constant  that  expresses  influence  of  temperature  on  sediment-water 
fluxes  (C?  l) 

T  =  temperature  (C °) 

b  =  constant  that  expresses  influence  of  light  on  sediment-water  fluxes 
(gm  m'2  langley1) 

Light  =  illumination  at  sediment-water  interface  Gangley  day1) 

By  convention  in  the  model,  fluxes  from  sediment  to  water  are  positive. 
Fluxes  from  water  to  sediment  are  negative.  The  empirical  model  indicates 
SOD  is  negative  (oxygen  flows  from  water  to  sediment)  when  no  light  is  avail¬ 
able  at  the  sediment-water  interface  (Figure  8-9).  The  magnitude  of  SOD 
increases  as  a  function  of  temperature.  The  availability  of  light  at  the  bottom 
makes  SOD  less  negative.  At  low  temperatures  and  high  light  intensity,  the 
bottom  gives  off  oxygen.  At  extremely  high  temperature,  however,  SOD  is 
negative  despite  light  availability.  (To  aid  in  interpretation  of  the  functions. 
Table  8-5  presents  a  few  reference  values  of  temperature  and  light) 
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SOD  (gm/m**2/day) 


Figure  8-9.  Effect  of  Temperature  and  Light  on  Sediment  Oxygen  Demand 


1  Table  8-5 

1 

Temperature  and  Light  at  Sediment-Water  Interface  at  Mean  Depth  of  Indian  I 
Rlver-Rehoboth  Bay  | 

■DTlul 

Tanparatura, 

Light,  Raglan  1.3, 
■npayoqr 

U0M,  Region  2,  1 
lanpayaay  B 

January 

53 

205 

46  . ”  | 

July 

26.6 

18 

*  1 

The  ammonium  flux  model  was: 


NH4flx  =  NH4flxb  ea(T"20)  -  bPN  Nlim  Light  (8-13) 


NH4flx  =  ammonium  flux  (gm  m'2  day*1) 

NH4flxb  =  ammonium  flux  at  20  C°  in  a  darkened  water  column 
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PN  =  algal  preference  for  ammonium 
Nlim  =  nitrogenous  nutrient  limitation  for  algae 
Definitions  of  remaining  terms  are  the  same  as  for  SOD  (Equation  8-12). 

Ammonium  flux  is  positive,  from  sediment  to  water,  when  no  light  is  avail¬ 
able  at  the  bottom  (Figure  8-10).  Sediment  ammonium  release  increases  as  an 
exponential  function  of  temperature.  When  light  is  available  at  the  bottom, 
ammonium  release  is  diminished.  Under  strong  illumination  at  low  tempera¬ 
ture,  sediments  take  up  ammonium.  The  nutrient  preference  and  limitation 
terms  are  identical  to  the  terms  specified  for  phytoplankton  (Equations  6-4  and 
6-21).  These  terms  prevent  sediment  ammonium  uptake  when  none  is  avail¬ 
able  in  the  water  column. 


FN  *  Nlim  *  Light  (langleys/day) 


Figure  8-10.  Effect  of  Temperature  and  Light  on  Sediment-Water  Ammonium  Flux 


Nitrate  moves  in  both  directions  across  the  sediment-water  interface.  This 
effect  is  influenced  by  concentration  in  the  water,  independent  of  benthic 
algae.  When  nitrate  is  abundant  in  the  water,  nitrate  flows  to  the  sediments 
where  it  is  denitrified  to  a  gaseous  form.  When  nitrate  is  absent  from  the 
water,  small  sediment  release  occurs  due  to  production  by  nitrification  in  the 
sediments  (Figure  8-11).  The  empirical  nitrate  model  was: 
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Figure  8-1 1 .  Effect  of  Temperature  and  Concentration  on  Sediment-Water  Nitrate  Flux 


N03flx  «  s(N03sed  -  N03wat)eaCr  '  20) 
-  (l  -  PN)Nlimb  Light 


(8-14) 


N03flx  =  nitrate  flux  (gm  m2  day'1) 

s  =  sediment-water  mass  transfer  coefficient  (m  day*1) 

N03sed  =  sediment  nitrate  concentration  at  sediment-water  interface 

(gm  m'3) 

N03wat »  water-column  nitrate  concentration  at  sediment-water  interface 
(gm  m'3) 

Definitions  of  remaining  terms  are  the  same  as  for  SOD  (Equation  8-12)  and 
ammonium  (Equation  8-13). 

No  influence  of  benthic  algae  on  nitrate  flux  could  be  discerned  in  the 
limited  data  base.  Potential  influence  of  benthic  algae  was  included  in  die 
empirical  model,  however.  The  influence  of  benthic  algae  on  nitrate  is 
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identical  to  the  influence  on  ammonium  except  that  ammonium  in  the  water 
must  be  first  exhausted  before  nitrate  is  consumed.  The  effect  of  benthic  algae 
is  to  increase  removal  of  nitrate  from  the  water  column  beyond  the  amount 
consumed  in  denitrification  (Figure  8-12).  The  enhanced  removal  occurs  only 
when  light  is  available,  however,  and  increases  in  proportion  to  light 
availability. 


Figure  8-12.  Effect  of  Concentration  and  Light  on  Sediment-Water  Nitrate  Flux 


The  phosphate  flux  model  was  identical  in  form  to  die  ammonium  flux 
model: 

P04flx  =  P04flxb  ea(T  _20)  -  bPlim  Light  (8-15) 

P04flx  =  phosphate  flux  (gm  m'2  day'1) 

P04flxb  =  phosphate  flux  at  20  C°  in  a  darkened  water  column 

Plim  =  phosphorous  nutrient  limitation  for  algae 

Definitions  of  remaining  terms  are  the  same  as  for  SOD  (Equation  8-12)  and 
ammonium  (Equation  8-13). 
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PCM  Flux  (mg/m**2/day) 


Phosphorus  flux  is  positive,  from  sedim,....  'o  water,  when  no  light  is  avail¬ 
able  at  the  bottom  (Figure  8-13).  Sediment  ^nosphate  release  increases  as  an 
exponential  function  of  temperature.  When  light  is  available  at  the  bottom, 
phosphate  release  is  diminished.  Under  strong  illumination  at  low  temperature, 
sediments  take  up  phosphate.  The  nutrient  limitation  term  is  identical  to  the 
term  specified  for  phytoplankton  (Equations  6-4).  This  term  prevents  sediment 
phosphate  uptake  when  none  is  available  in  the  water  column. 


Plim  *  Light  (langleys/day) 


Figure  8-13.  Effect  of  Temperature  and  Light  on  Sediment-Water  Phosphate  Flux 


Parameters  in  the  empirical  models  were  initially  evaluated  through  nonlin¬ 
ear  regression.  Initial  values  were  refined  for  uniformity  and  to  optimize 
model  calibration.  Final  values  are  presented  in  Table  8-6.  Sediment-water 
fluxes  determined  by  the  functions  are  compared  to  observations  in 
Figure  8-14. 


Computation  of  BOD 

The  biochemical  oxygen  demand  (BOD)  analysis  combines  measures  of  the 
concentration  of  oxygen-demanding  material  and  the  oxidation  rate.  BOD  is 
not  a  model  state  variable  but  the  BOD  of  quantities  computed  by  the  model 
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Table  M 

Parameters  In  Empirical  Flux  Models 


Baas  Flux, 
gm  m2  day'1 

B 

b, 

gm  m*  langley'1 

mdey'1 

£0 

Qm  If  III 

Sediment 

Oxygen 

Demand 

-1.45 

0.05 

0.011 

Ammonium 

0.07 

0.05 

4.5X10"4 

Nitrate 

0.05 

4JX10-4 

0.1 

0.03 

Phosphate 

0.009 

0.05 

l.lxicr4 

Figure  8*14.  Observed  and  Modeled  Sediment-Water  Fluxes 

Chapter  8  Additional  Model  Inputs  and  Outputs 


1 

I 


can  be  calculated  for  comparison  to  five-day  BOD  observations.  Model  BOD3 
was  computed: 


BOD,  *  AOCRB (1  -  e‘  5BM)  ♦  AOCRDOC(l  -  e'  5Kdoc) 

(8-16) 

-  5  NTm 

♦  AONTNH4(l  -  e  KKnnt) 


BODs  =  five-day  biochemical  oxygen  demand  (gm  m'3) 

A  OCR  =  dissolved  oxygen-to- carbon  ratio  in  respiration 
(2.67  gm  C>2  gm  C1) 

B  =  algal  biomass  (gm  C  m3) 

BM  =  basal  metabolism  (day1) 

DOC  =  dissolved  organic  carbon  (gm  m'3) 

Kdoc  =  respiration  rate  of  dissolved  organic  carbon  (day1) 

AONT  =  mass  dissolved  oxygen  consumed  per  mass  ammonium  nitrogen 
nitrified  (4.33  gm  Oj  gm'1  N) 

NTm  =  maximum  nitrification  rate  (gm  N  m'3  day'1) 

KHnnt  =  half-saturation  concentration  of  NH4  required  for  nitrification 
(gm  N  m'3) 

The  relationship  is  an  adaptation  of  die  formulae  for  oxygen  consumption 
expressed  in  Chapter  6.  In  the  computation  of  BOD,  oxygen  is  considered  to 
be  available  in  excess  and  ammonium  is  considered  to  be  much  less  than  the 
half-saturation  concentration,  KHnnt.  Since  the  BOD  analysis  is  conducted  in 
darkness,  algal  respiration  is  considered  but  not  photosynthesis. 


Diurnal  Dissolved  Oxygen  Fluctuations 

Algal  physiological  functions  that  create  diurnal  dissolved  oxygen  fluctua¬ 
tions  are  not  well-represented  by  the  conventional  photosynthesis  and  respira¬ 
tion  relations  employed  in  eutrophication  modeling.  Additional  influences  on 
diurnal  fluctuations,  particularly  macrophytes,  are  not  represented  at  all  in  the 
model.  Consequently,  the  preferred  means  of  computing  daily-minimum  dis¬ 
solved  oxygen  is  through  a  relationship  to  predicted  daily-mean  dissolved 
oxygen.  The  relationship  was  developed  by  fitting  a  sinusoidal  function  (Fig¬ 
ure  8-13)  to  each  diurnal  dissolved  oxygen  measure: 
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Figure  8-15.  Sinusoidal  Fit  to  Diurnal  Dissolved  Oxygen  Data 


DO 


*CC°Y5T 


2n 


(8-17) 


DO  =  observed  dissolved  oxygen  concentration  (gm  02  m'3) 

a  =  mean  dissolved  oxygen  concentration  (gm  02  m'3) 

b,  c  =  parameters  that  determine  amplitude  and  phase  of  diurnal  DO 
fluctuation  (gm  02  m"3) 

Parameters  a,  b,  and  c  were  determined  by  regression.  Amplitude  was 
computed: 


A  *  /b2  ♦  c2 


(8-18) 


A  =  amplitude  of  diumal  DO  fluctuation  (gm  02  m'3) 
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Next,  regression  was  used  to  relate  amplitude  to  chlorophyll  concentration. 
The  resulting  relationship  was: 


A  =  1.51  ♦  0.033  Chi 


(8-19) 


Chi  =  chlorophyll  ’a’  concentration  (mg  m'3) 


R2  for  the  regression  was  0.76.  Worth  noting  is  that  a  diurnal  fluctuation  of 
1.5  gm  Oj  m'3  exists  in  the  absence  of  algae  (Figure  8-16)  and  is  presumably 
caused  by  aquatic  vegetation. 

To  compute  daily-minimum  DO,  daily-average  predicted  chlorophyll  is 
substituted  in  Equation  8-19  to  obtain  amplitude.  Amplitude  is  subtracted 
from  daily-average  predicted  dissolved  oxygen  to  yield  daily  minimum. 


Figure  8-16.  Effect  of  Chlorophyll  on  Amplitude  of  Diurnal  Dissolved  Oxygen 
Fluctuation 
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Chapter  IX:  Water  Quality 
Model  Application  and 
Analysis 


Water  Quality  Time  Series 

The  primary  means  of  calibrating  the  water  quality  model  was  through 
comparison  of  modeled  and  observed  time  series  of  water  quality  constituents. 
Time  series  were  compared  in  eight  model  cells  (Figure  9-1)  selected  to  pro¬ 
vide  spatial  coverage  and  a  sufficient  number  of  observations.  Calibration  was 
an  iterative  process  in  which  kinetics  coefficients  (Table  6-2),  sediment  model 
parameters  CTable  8-6),  and  load  fractions  (Table  8-2)  were  adjusted  to 
improve  model-data  agreement  Calibration  ceased  when  no  additional 
improvement  was  possible.  The  final  calibration  time  series  are  shown  as 
Figures  9-2  through  9-9.  Observations  are  grab  samples  collected  within  a 
model  cell.  Symbols  indicate  the  source  of  the  observations.  Model  results 
are  daily  average  values  for  the  cell.  The  temporal  origin  (Year  0)  is  Janu¬ 
ary  1,  1988.  Time  series  cease  (Year  3)  on  December  31,  1990. 


Longitudinal  Transect 

Additional  insight  into  the  calibration  status  was  gained  through  comparison 
of  observations  and  model  results  along  a  longitudinal  transect  that  extended 
from  the  Millsboro  spillway  to  the  mouth  of  the  Indian  River  inlet  Transects 
were  plotted  for  the  summer  of  each  calibration  year  (Figures  9-10  through 
9-12).  The  spatial  origin  (Km.  0)  was  at  the  spillway.  Data  collection  along 
the  transect  was  not  synoptic.  Neither  were  observations  collected  at  uniform 
intervals  throughout  the  summer.  Consequently,  no  attempt  was  made  to  pro¬ 
duce  a  summer-average  from  the  data.  Grab  samples  collected  in  die  months 
of  June,  July,  or  August  were  plotted  at  the  location  of  the  observation.  These 
were  compared  to  summer-average  model  results.  Point-by-point  agreement 
between  the  model  average  and  grab  samples  was  not  expected.  Rather,  die 
{dots  were  used  to  assess  observed  and  predicted  spatial  trends  and  approxi¬ 
mate  agreement  in  concentration. 
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Figure  9-1 .  Grid  Cells  Selected  for  Time  Series  Presentation 


Sediment-Water  Fluxes 

Time  series  of  modeled  sediment- water  fluxes  (Figures  9-13  through  9-20) 
were  produced  at  the  same  locations  as  the  calibration  time  series  (Figure  9-1). 
By  convention,  fluxes  from  water  to  sediments  were  shown  as  negative  values. 
Positive  values  indicated  sediment  release  to  the  water  column.  The  “Carbon  - 
Oxygen  Flux”  time  series  included  carbon  deposition  to  the  sediments  and 
sediment  oxygen  demand.  The  “Nitrogen  Flux”  time  series  included  particu¬ 
late  nitrogen  deposition  to  the  sediments,  sediment-water  ammonium  flux,  and 
sediment-water  nitrate  flux.  The  “Phosphorus  Flux”  time  series  indicated 
particulate  phosphorus  deposition  to  the  sediments  and  sediment-water  phos¬ 
phate  flux. 


Interpretive  Information 

Time  series  of  interpretive  information  available  from  the  model  (Figures  9- 
21  through  9-28)  were  produced  at  the  same  locations  as  the  calibration  lime 
series  (Figure  9-1).  One  time  series  indicated  the  computed  nutrient  and  light 
limitations  on  algal  growth  (Equations  6-4,  6-6).  Potential  limitation  ranged 
from  zero  (complete  limitation  on  growth)  to  unity  (no  limit  on  growth).  A 
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Series  of  Observed  and  Modeled  Water-Column  Concentrations,  Grid  Cell  158  (SI 
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Figure  9-2.  (Sheet  2  of  3) 


9-6 


Chapter  9  Water  Quality  Model  Application  and  Analysis 


Figure  9-3.  Time  Series  of  Observed  and  Modeled  Water-Column  Concentrations,  Grid  CeH  163  (Sheet  1  ot  3) 
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Figure  9*4.  (Sheet  2  of  3) 
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Figure  9-5.  Time  Series  of  Observed  and  Modeled  Water-Column  Concentrations,  Grid  Cell  143  (Sheet  1  of  3) 
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Figure  9-5.  (Sheet  3  of  3) 
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Figure  9-6.  Time  Series  of  Observed  and  Modeled  Water-Column  Concentrations,  Grid  CeR  203  (Sheet  1  of  3) 
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Figure  9-6.  (Sheet  2  of  3) 
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Figure  9-7.  Time  Series  of  Observed  and  Modeled  Water-Column  Concentrations,  Grid  Cell  220  (Sheet  1  of  3) 


Observed  and  Modeled  Water-Column  Concentrations, 
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Figure  9-8.  (Sheet  2  of  3) 
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Figure  9-9.  Time  Series  of  Observed  and  Modeled  Water-Column  Concentrations,  Grid  Cell  275  (Sheet  1  of  3) 
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Figure  9-9.  (Sheet  3  of  3) 
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Figure  9-10.  Observed  and  Modeled  Water-Column  Concentrations  Along  Transect  From  Millsboro  to  Indian  River  Inlet,  Summer 
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Figure  9-1 1 .  (Sheet  3  of  3) 


Figure  9-12.  (Sheet  2  of  3) 
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Figure  9-15.  Time  Series  of  Modeled  Sediment- Water  Fluxes,  Grid  Cell  116 
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Figure  9*17.  Time  Series  of  Modeled  Sediment-Water  Fluxes,  Grid  Cell  203 


Chapter  9  Water  Quality  Modal  Application  and  Analysis 


Looi 


L%' 

^  Ji  o 

L'-J 

ccb-S: 


SR°R^R8 


tn 


ui 

>- 


A0/2**W/9W 


se 

CM 

© 

o 


o 


1 

I 


I 


O 

9 


9-43 


9-44 


Chapters  Water Quality Modal Appfcation and Anatytia 


Figure  9-21.  Time  Series  of  Interpretive  Information,  Grid  Cell  158 


Series  of  Interpretive  Information,  Grid  Cel  163 
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Figure  9-23.  Time  Series  of  Interpretive  Information,  Grid  Cel  1 16 
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Figure  9-25.  Time  Series  of  Interpretive  Information,  Grid  Celt  203 
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Figure  9-27.  Time  Series  of  Interpretive  Information,  Grid  Cell  252 


second  plot  indicated  net  areal  primary  production,  computed  as  algal  carbon 
production  minus  total  water  column  carbon  consumption.  A  third  plot  indi¬ 
cated  instantaneous  water-column  respiration,  on  a  volumetric  basis.  Computa¬ 
tion  of  respiration  was  analogous  to  computation  of  BOD  (Equation  8-16) 
except  no  five-day  integration  was  performed.  Together,  the  primary  produc¬ 
tion  and  respiration  plots  indicated  biological  activity.  Regions  with  high 
production  and  respiration  are  more  eutrophic  than  regions  with  lesser  produc¬ 
tion  and  respiration.  The  final  time  series  was  modeled  light  extinction,  the 
sum  of  background  light  extinction  input  to  the  model  and  attenuation  due  to 
computed  chlorophyll  concentration. 


Evaluation  of  Model  Performance 

Graphical  Evaluation 

Evaluation  of  model  performance  from  a  large  number  of  time  series  is 
difficult.  To  aid  in  evaluation,  graphical  and  statistical  summaries  were  per¬ 
formed.  For  major  water  quality  constituents,  scatterplots  of  observed  versus 
modeled  concentrations  were  produced  (Figure  9-29).  The  summary  included 
each  observation  in  the  data  base,  not  only  observations  in  the  time  series.  As 
with  the  time  series,  grab  samples  were  plotted  against  daily  averages  in  model 
cells.  A  solid  line  indicating  one-to-one  correspondence  was  shown  on  each 
plot.  Cumulative-difference  plots  were  also  produced.  These  indicated  the 
percentage  of  observations  which  differed  from  modeled  concentrations  by  less 
than  an  indicated  value. 

The  temperature  scatterplot  indicated  ideal  model-data  agreement.  Discrep¬ 
ancies  between  predictions  and  observations  were  small.  Observations  and 
modeled  temperatures  fell  nearly  on  a  one-to-one  line  throughout  the  observed 
range. 

Salinity  correspondence  was  also  excellent,  especially  in  the  upper  salinity 
range.  Below  =20  ppt,  the  model  showed  a  tendency  to  underestimate  salinity, 
however. 

The  chlorophyll  scatterplot  was  typical  for  this  substance.  Observations  and 
model  agreed  well  in  lower  chlorophyll  ranges.  Extreme  observations,  above 
SO  mg  m'\  existed  which  could  not  be  modelled,  however.  In  Indian  River/ 
Rehobo th  Bay,  most  of  these  observations  were  collected  near  Millsboro.  The 
observations  were  near-surface  samples  and  resulted  from  algal  mat  formation 
and  concentration  of  the  mats  by  wind  and  physical  features.  The  model, 
which  computes  depth-average  concentrations,  is  not  suited  for  prediction  of 
chlorophyll  on  the  scale  of  these  observations. 

The  BOD  scatterplot  indicated  correct  model  representation  of  the  processes 
which  produce  the  range  of  observed  BOD.  Observed  and  modeled  maxima 
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Figure  9-29.  Graphical  and  Statistical  Comparison  of  Individual  Observations  Versus  Model 
Results  (Sheet  1  of  6) 


Chapter  9  Water  Quality  Model  Application  and  Analysis 


9-53 


observed 


Chlorophyll 


y  lntere*pl-26.388 
r««2  -  0.137 
•lop*  ■  3.271 


6005 

y  Intercept  1 .884 
r«w2  -  0.477 
•lop*  -  0.862 


o 


mot  *rror  -  15.897 
absolute  ■•an  *rror  -  30.925 
RMS  -  102.257 
r*lotlw*  *rror  *  89.62 


abs ( observed-predicted ) 


Man  error  -  1.406 
absolute  Man  error  •  1.521 
RMS  -  1.852 
relative  error  -  31.2% 


obs ( observed-predicted ) 


Figure  9-29.  (Sheet  2  of  6) 


9-54 


Chapter#  Water  Quteity  Model  Appicabon  and  Analysis 


observed  observed 

)  2.0  3.0  4.0  0.00  0.25  0.r0  0.75  1.00 


Rinmontum 


y  Intercept  0.078 
p«h 2  ■  0.032 
•lop*  "  0.173 


Nitrate-Nitrite 


y  Intercept  0.064 
re* 2  -  0.727 
elope  -  0.687 


Figure  9*29.  (Sheet  3  of  6) 


neon  error  ”"0.016 
absolute  mean  error  -  0.094 
RMS  -  0.145 

relative  error  ■  96.82 


0.00  0.25  0.50  0.75 

abs ( observed-predicted ) 


•eon  error  - -0.037 
abeolute  neon  error  -  0.173 
RMS  -  0.338 
relative  error  ■  49.62 


0  12  3 

abs ( obser ved-predlcted ) 


Chapter  9  Water  Quality  Model  Application  and  Analysis 


9-55 


Total  Kjeldahl  N 

y  intercept  0.072 
r«*2  -  0.365 
•lop*  -  0.983 


■won  error  -  0.055 
absolute  naan  «rror  -  0.332 
RMS  -  0.453 
relative  •rror  “  31. 4  X 


1.0  2.0  3.0 

predicted 


0  0.5  1  1.5  2 

abc ( obs«rv®d_predlcted ) 


Total  Nitrogen 

jj  Intercept  0.297 
r**2  -  0.558 
•Lop*  -  0.790 


■wan  •rror  -  0.030 
absolute  wean  error  -  0.357 
RrtS  -  0.487 

relative  error  -  27.5 X 


0.0  2.5  5.0 

predicted 


Figure  9-29.  (Sheet  4  of  6) 


abs ( observed-predicted ) 


Chapters  Water  Quality  Model  AppScaion  and  Analysis 


Dissolved  Inorganic  P 


y  intorcspl  0.022 
r«N2  -  0.001 
slops  “  0.067 


Total  Phosphorus 

y  Intsrcspt  0.023 
p*m2  -  0.169 
slops  “  0.944 

«o 


•son  srror  -  0.010 
absoluts  n>s an  srror  ■  0.021 
RMS  -  0.048 
rslotlvs  srror  “  91.3% 


abs ( obser ved-pr®dlc t#d ) 


•son  srror  •  0.019 
absoluts  nsan  srror  *  0.035 
RMS  -  0.061 

rstotlvs  srror  “  36.4% 


abs ( obser ved~predlct8d ) 


Figure  9-29.  (Sheet  5  of  6) 


Chapter  9  Water  Quality  Modal  Appicadon  and  Analysis 


5.0  10.0  15.0 


0c«*olv«d  Oxygon 

y  Intercept  3.652 
r""2  ~  0.024 
•lop*  ■  0.283 


Mon  trror  “”1.411 
otiolul*  Man  *rror  •  2.079 
RMS  -  2.750 
r*lotLV*  *rror  -  36.8% 


aba ( obeerved~predLcted ) 


Figure  9-29.  (Sheet  6  of  6) 


9-58 


Chapters  Water Qualfy Modal AppficaSon and Analysis 


and  minima  coincided.  Modeled  BOD  was  usually  less  than  observed, 
however. 

The  ammonium  scatterplot  was  typical  for  this  substance.  The  range  of 
modeled  ammonium  agreed  with  the  observed  range  but  little  one-to-one  corre¬ 
spondence  between  model  and  observations  existed.  The  nitrate  scatterplot 
indicated  good  correspondence  between  model  and  observations  throughout  the 
observed  range.  Good  model-data  correspondence  also  existed  for  total 
Kjeldahl  nitrogen  and  total  nitrogen  although  extreme  total  Kjeldahl  measures 
were  not  replicated  by  the  model. 

The  phosphate  scatterplot  was  similar  to  ammonium.  The  range  of  modeled 
phosphate  agreed  with  the  observed  range  but  little  one-to-one  correspondence 
existed,  especially  for  the  extreme  values.  Correspondence  between  modeled 
and  observed  total  phosphorus  was  better  than  for  phosphate.  Only  a  few 
extreme  observations  were  not  replicated. 

The  dissolved  oxygen  scatterplot  indicated  behavior  similar  to  ammonium 
and  phosphate.  Observed  and  modeled  range  were  reasonable  but  point-by¬ 
point  comparisons  showed  large  discrepancies.  The  lowest  observations  were 
not  well-replicated.  One  reason  for  the  discrepancies  is  the  comparison  of 
daily-average  model  values  with  instantaneous  observations.  Since  the  obser¬ 
vations  are  subject  to  diurnal  fluctuations,  direct  correspondence  cannot  be 
expected.  The  minimum  values  in  the  observations  are  better  represented  by 
application  of  the  diurnal  curve  technique  described  in  Chapter  8  than  by  direct 
utilization  of  model  output. 

The  large  number  of  observations  partially  obscures  interpretation  of  the 
scatterplots.  The  comparison  of  individual  observations  also  implies  temporal 
and  spatial  resolution  that  are  not  available  in  the  model.  The  grab  samples 
are  point  observations  influenced  by  physical  processes  that  occur  on  time 
scales  less  than  a  tidal  cycle.  The  smallest  spatial  scale  resolved  by  the  model 
is  the  cell  length,  roughly  500  m.  While  the  model  time  step,  roughly  five 
minutes,  is  consistent  with  the  temporal  scale  of  processes  that  influence  the 
observations,  forcing  functions  such  as  loads  and  light  extinction  are  known 
only  on  a  monthly  basis.  Improved,  mote  reliable  views  of  model  performance 
are  obtained  by  temporal  and  spatial  aggregation  of  observations  and  model 
output  so  that  small-scale,  essentially  random,  variations  are  removed  from  the 
observation  set 

Aggregation  was  performed  by  regional  averaging  of  all  data  collected  in 
the  months  of  June,  July,  and  August.  These  months  were  selected  since  they 
contained  most  data  and  since  summer  is  the  season  of  critical  water  quality. 
Regions  were  the  same  as  defined  for  analysis  of  light  extinction  (Figure  8-7). 
An  average  was  computed  for  each  region  for  summer  of  1988, 1989,  and 
1990.  (For  some  season-region  combinations,  no  observations  existed,  how¬ 
ever.)  The  averaged  observations  were  compared  to  model  results  averaged 
over  identical  regions  and  seasons  (Figure  9-30).  The  analysis  also  permitted 
evaluation  of  model  performance  on  a  regional  basis. 
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Figure  9-30.  Graphical  and  Statistical  Comparison  of  Summer-Average  Observations  Versus 
Model  Results  (Sheet  1  of  6) 
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Modeled  temperature  was  ideal  on  a  point- by-point  basis  so  little  additional 
information  was  gained  by  examination  of  aggregated  values.  Aggregation  of 
salinity,  however,  indicated  nearly  perfect  model-data  agreement  for  six  of  nine 
seasonal-regional  averages.  Model  chlorophyll  results  were  perfect  in  Indian 
River  Bay.  The  model  overestimated  chlorophyll  in  Rehoboth  Bay,  however. 
Chlorophyll  in  upper  Indian  River  was  underestimated  for  reasons  previously 
explained.  The  behavior  of  aggregated  BOD  was  the  same  as  for  individual 
comparisons.  The  model  indicated  the  correct  trend  in  BOD  but  underesti¬ 
mated  concentration.  Large  portions  of  the  erratic  ammonium  comparisons 
were  removed  by  seasonal-regional  averaging.  Model-data  ammonium  com¬ 
parisons  were  ideal  for  five  of  nine  seasonal-regional  combinations.  Erratic 
predictions  and  observations  still  existed  for  some  seasons  and  regions,  how¬ 
ever.  Nitrate,  total  Kjeldahl  nitrogen,  and  total  nitrogen  comparisons  on  the 
seasonal-regional  basis  mirrored  the  favorable  performance  of  individual  com¬ 
parisons.  Phosphate  demonstrated  a  clear  trend  of  underprediction  by  the 
model  in  all  regions.  Total  phosphorus  demonstrated  ideal  agreement  of  pre¬ 
dictions  and  observations  in  most  regions  and  seasons  although  the  model 
underestimated  in  some  comparisons.  Comparisons  indicated  the  model  over¬ 
estimated  dissolved  oxygen  in  most  season-region  combinations  but  correspon¬ 
dence  in  several  was  near-perfect. 


Statistical  Evaluation 

The  scatterplots  are  useful,  qualitative  summaries  of  model  performance. 
Quantitative  measures  are  preferred  by  some  analysts,  however.  No  consensus 
exists  on  a  standard  suite  of  quantitative  performance  measures.  As  a  result, 
we  compute  a  large  suite,  reported  on  the  same  figures  as  the  scatterplots 
(Figures  9-29,  9-30). 

Regression  Analysis  -  One  popular  form  of  model-data  comparison  treats 
model  results  as  independent  variables  and  observations  as  dependent  variables 
in  linear  regression  analysis.  A  high  degree  of  correspondence  between  model 
and  observations  is  indicated  by  a  correlation  coefficient  near  unity,  an  inter¬ 
cept  near  zero,  and  a  slope  near  unity. 

Mean  Error  -  The  mean  error  statistic  is  defined: 


ME  = 


£(0  -P) 

n 


(9-1) 


ME  =  mean  error 
O  =  observation 
P  =  model  result 
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n  =  number  of  observations 


The  mean  error  is  a  summary  of  the  model  tendency  to  overestimate  or 
underestimate  the  data.  The  mean  error  can  be  zero  even  though  large  discrep¬ 
ancies  exist  in  individual  model-data  comparisons.  The  mean  error  is  espe¬ 
cially  useful  to  indicate  potential  errors  in  loads,  sediment-water  fluxes,  and 
other  sources  and  sinks  to  the  water  column. 

Absolute  Mean  Error  -  The  absolute  mean  error  statistic  is  defined: 


AME 


Ei°  -pl 

n 


(9-2) 


AME  =  absolute  mean  error 

The  absolute  mean  error  is  a  measure  of  the  average  discrepancy  between 
observations  and  model  results.  The  absolute  mean  error  does  not  differentiate 
between  overestimates  and  underestimates  by  the  model. 

Root  Mean  Square  Error  -  The  root  mean  square  error  statistic  is  defined: 


RMS 


E<°  -P)2 


N 


n 


(9-3) 


RMS  =  root  mean  square  error 

The  root  mean  square  error  is  an  alternate  quantification  of  the  average 
discrepancy  between  observations  and  model  results. 

Relative  Error  -  Tire  relative  error  is  defined: 


RE  = 


El<)  -pi 


(9-4) 


RE  =  relative  error 

The  relative  error  statistic  normalizes  absolute  mean  error  by  the  magnitude 
of  the  observations.  The  relative  error  statistic  is  especially  useful  in  compar¬ 
ing  performance  between  variables  of  different  magnitude  or  performance  of 
models  of  different  systems. 
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The  mean  error  and  root  mean  square  error  for  major  water  quality  constit¬ 
uents  are  summarized  in  Table  9-1.  Summaries  are  for  instantaneous  and 
seasonal-regional  comparisons. 


Table  9-1 

Summary  of  Mean  and  Root  Mean  Square  Error 

mean  error 

mssn  error 

DM6  eifinr 
IDD  tnOf 

RMS  error  [ 

■ - ■  —  fillies 

rasn  Knfif| 
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\Mw  POtfuS 

inflWn  KIW| 

Summer  Regional 
Average 

■ — ■ —  f|i, 

■osn  HirPi 

Individual 
Date  Points 

-  I 

■KMn  DDf)  1 

Summer  Regional 
Average  j 

Temperature, 

C° 

-0.08 

0.84 

225 

1.49 

Salinity, 

PPt 

2.53 

-2.16 

5.9 

3.46 

Chlorophyll, 
mg  nr* 

15.9 

29.1 

102.3 

74.1 

BOOS, 
gm  m'5 

Ml 

237 

1.85 

2.37 

Ammonium, 
gm  m~5 

-0.016 

0.029 

0.145 

0.115 

Nitrate, 
gm  m~3 

-0.037 

0.102 

0.338 

0.111 

Total  Kjeldahl 
N,  gm  m~J 

0.055 

0.112 

0.453 

0283 

Total  N, 
gm  m~J 

0.03 

0209 

0.487 

0.344 

Phosphate, 

gm  nrf5 

0.01 

0.021 

0.048 

0.022 

Total  P, 
gm  m"1 

0.019 

0.024 

0.061 

0.035 

Dissolved 
Oxygen, 
gm  m"J 

-Ml 

-122 

2.75 

1.77 

The  summary  indicates  mean  error  is  generally  larger  in  magnitude  for  the 
summer  average  than  for  individual  observations.  One  interpretation  of  this 
phenomenon  is  that  sources  and  sinks  for  summer  only  are  quantified  to  a 
lesser  degree  of  accuracy  than  for  the  time  period  encompassed  by  the  entire 
data  base  (April  through  September).  The  RMS  error  for  the  seasonal-regional 
average  is  about  half  the  error  for  instantaneous  comparisons.  The  implication 
of  this  statistic  is  that  the  model  is  more  accurate  in  reproducing  spatial  and 
temporal  aggregations  than  individual  observations. 
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Interpretation  of  performance  statistics  is  facilitated  when  some  basis  for 
comparison  is  present.  We  compared  the  relative  error  statistics  with  results  of 
the  Chesapeake  Bay  model  study  (Cerco  and  Cole  1993).  Individual  and  sea¬ 
sonal-regional  averages  from  Indian  River  were  compared  to  seasonal-regional 
averages  from  the  Bay  (Figure  9-31).  The  Chesapeake  Bay  study  made  use  of 
an  extensive,  high-quality  data  base  and  employed  highly  accurate  measures  of 
loads  and  boundary  conditions.  Calibration  of  the  model  consumed  several 
years  and  was  subject  to  extensive  review.  We  regard  the  performance  statis¬ 
tics  of  the  Bay  study  as  current  state  of  the  art  for  long-term  applications  of 
comprehensive  eutrophication  models. 

Relative  error  for  seasonal-regional  averages  in  the  Indian  River  model 
application  was  generally  greater  than  in  Chesapeake  Bay  but  comparable  for 
all  but  two  substances.  Relative  error  for  chlorophyll  in  Indian  River  was 
double  the  Chesapeake  Bay  value.  The  chlorophyll  statistic  in  Indian  River 
reflected  the  extreme  concentrations  observed  at  the  surface  near  the  Millsboro 
spillway.  Concentrations  in  excess  of  1000  mg  m'3  were  in  the  data  base. 
Discrepancies  between  these  observations  and  model  results  added  enormous 
bias  to  the  performance  statistics.  Relative  error  for  phosphate  was  also  double 
the  Chesapeake  Bay  value.  We  were  informed  by  DNREC  that  the  phosphate 
data  in  STORET  is  dubious.  We  noted  the  STORET  phosphate  data  fre¬ 
quently  disagreed  with  the  observations  collected  by  the  College  of  Marine 
Studies  (Figures  9-3,  9-4,  9-7,  9-8).  Most  often,  when  conflicts  in  the  data 
were  noted,  model  results  conformed  to  die  College  of  Marine  Studies  rather 
than  the  STORET  data.  We  concluded  the  large  phosphate  relative  error  was 
due  to  faulty  data  rather  than  a  shortcoming  in  the  model  calibration. 


The  Annual  Cycle  of  Water  Quality 

The  time  series  of  conditions  in  the  water  column,  sediment-water  interac¬ 
tions,  and  interpretive  quantities,  together  with  other  information  available 
from  the  model  permit  a  description  of  the  annual  cycle  of  water  quality  in 
Indian  River/Rehoboth  Bay. 

During  the  first  four  months  of  the  year,  chlorophyll  is  virtually  absent 
from  die  water  column  (Figures  9-2  through  9-9).  Nitrogen,  in  die  form  of 
nitrate,  and  light  to  support  algal  growth  are  readily  available  (Figures  9-21 
through  9-28).  A  slight  phosphorus  limitation  to  growth  is  present  but  the 
primary  influence  on  algal  abundance  is  temperature.  The  water  is  too  cold  to 
support  substantial  primary  production  by  phytoplankton.  In  May,  the  water 
warms  sufficiendy  to  support  algal  production  so  that  chlorophyll  concentra¬ 
tion  rapidly  rises  to  its  summer  maximum  value.  Limitations  to  algal  growth 
during  summer  vary.  Light  extinction,  due  to  background  extinction  and  algal- 
self  shading,  is  at  its  annual  peak  (Figures  9-21  through  9-28).  In  upper 
Indian  River,  phosphorus  is  die  limiting  nutrient  (Figures  9-21  through  9-23). 
In  Lower  Indian  River  and  in  Rehoboth  Bay,  however,  nitrogen  is  the  domi¬ 
nant  nutrient  limit  on  algal  growth  (Figures  9-24  through  9-28).  Chlorophyll 
concentration  fluctuates,  in  response  to  flow  and  loading  events,  until 
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Figure  9-31.  Relative  Error  in  Indian  RKrer/Rehoboth  Bay  and  Chesapeake  Bay  Model  Applications 


mid-October.  In  October,  temperature  falls  below  the  level  needed  to  support 
substantial  algal  growth  and  chlorophyll  declines  to  its  winter  minimum 
concentration. 

During  the  first  four  months  of  the  year,  background  light  extinction  is  near 
its  minimum  value  and  light  attenuation  by  chlorophyll  is  negligible  (Fig¬ 
ures  9-21  through  9-28).  Light  is  available  to  support  production  by  benthic 
algae.  As  a  consequence,  dissolved  oxygen  is  produced  at  the  sediment-water 
interface  and  released  to  the  water  column  (Figures  9-13  through  9-20). 
Ammonium  and  phosphate  are  removed  from  the  water  column  and  consumed 
at  die  sediment-water  interface.  Nitrate  transport  is  into  the  sediments,  driven 
by  the  high  concentrations  in  the  overlying  water.  Two  phenomenon  accom¬ 
pany  the  wanning  of  the  water  column  in  May.  First,  oxygen  consumption 
and  nutrient  production  in  the  subsurface  sediments  increase.  Second,  light 
extinction  increases  due  to  the  seasonal  variation  in  background  extinction  and 
chlorophyll  attenuation.  As  a  consequence  of  temperature-enhanced  sediment 
processes  and  shading  of  the  benthic  algae,  the  sediments  consume  oxygen 
removed  from  the  water  column.  The  sediments  become  sources  of  ammon¬ 
ium  and  phosphate  to  the  water.  Sediment-water  nitrate  flux  is  negligible  due 
to  diminished  concentration  in  the  overlying  water.  Sediment  consumption  of 
oxygen  and  production  of  ammonium  and  phosphate  continue  until  autumn 
cooling  of  the  water  column.  Cooling  diminishes  diagenetic  processes  in  the 
sediments  and  is  accompanied  by  diminished  light  extinction.  Benthic  algae 
become  active  so  that,  by  the  end  of  the  year,  dissolved  oxygen  is  again  pro¬ 
duced  at  the  sediment-water  interface.  Ammonium  and  phosphate  enter  a 
transition  from  sediment  release  to  consumption  at  the  sediment-water  inter¬ 
face. 

The  annual  nitrogen  cycle  varies  according  to  locatioa  In  lower  Indian 
River  and  Rehoboth  Bay  (Figures  9-5  through  9-9),  total  nitrogen  concentra¬ 
tion  demonstrates  a  sinusoidal  behavior  with  maximum  concentration  in  late 
summer.  The  phase  of  the  total  nitrogen  function  corresponds  to  the  phase  of 
annual  temperature  variation.  The  correspondence  occurs  because  temperature 
drives  sediment  ammonium  release  (Figures  9-13  through  9-20),  a  primary 
source  of  nitrogen  to  the  lower  Bays.  In  upper  Indian  River,  total  nitrogen  is 
more  closely  related  to  runoff  events  so  that  maximum  concentrations  typically 
occur  in  late  winter  or  early  spring  (Figures  9-2  through  9-4).  Since  nitrate  is 
a  major  constituent  of  runoff  nitrogen,  nitrate  maxima  occur  in  winter  and 
spring  throughout  the  system.  The  magnitude  of  the  maximum  diminishes 
away  from  the  runoff  entry  points,  however.  In  the  lower  Bays,  in  summer, 
nitrate  is  virtually  absent  due  to  algal  uptake  and  negligible  local  runoff.  The 
ammonium  cycle  is  complex,  influenced  by  external  loading,  sediment 
production/consumption,  and  algal  uptake.  In  much  of  the  system,  a  summer 
maximum,  coincident  with  sediment  release  would  be  expected.  As  rapidly  as 
ammonium  is  released  from  the  sediment,  however,  it  is  consumed  by  algae. 

In  lower  Indian  River  and  Rehoboth  Bay  (Figures  9-5  through  9-9)  maximum 
ammonium  concentrations  are  evident  in  October  through  December.  In  these 
months,  phytoplankton  uptake  is  diminished  but  uptake  by  benthic  algae  is  not 
sufficient  to  remove  the  ammonium  from  the  water  column.  Upper  Indian 
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River  also  demonstrates  die  autumn  peak  but  high  concentrations,  due  to  sedi¬ 
ment  release,  occur  in  summer  as  well.  In  this  portion  of  the  system,  summer 
algae  populations  tend  to  be  limited  by  phosphorus  availability  so  that  not  all 
ammonium  is  removed  by  algal  activity. 

The  annual  cycle  of  total  phosphorus  also  coincides  with  the  temperature 
cycle.  As  with  nitrogen,  the  coincidence  is  due  to  temperature-enhanced  nutri¬ 
ent  release  from  the  sediment.  The  influence  of  runoff  in  the  annual  phospho¬ 
rus  cycle  is  not  as  evident  as  for  nitrogen  however  (Figures  9-2  through  9-4). 

Annual  cycling  in  phosphate  is  similar  to  ammonium.  Lower  Indian  River  and 
Rehoboth  Bay  demonstrate  an  autumn  peak  (Figures  9-5  through  9-9)  that 
presumably  is  due  to  diminished  phytoplankton  uptake  coupled  with  minimal 
activity  by  benthic  algae.  Near  the  Mills  boro  spillway,  however,  phosphate 
exhibits  erratic  behavior  due  to  loading  from  runoff  events,  consumption  and 
release  in  sediments,  and  algal  uptake. 

As  with  most  estuarine  systems,  minimum  dissolved  oxygen  in  Indian 
River/Rehoboth  Bay  occurs  during  late  summer  when  water  temperature  is 
highest.  The  coincidence  of  high  temperature  and  low  dissolved  oxygen  is  due 
to  temperature-enhanced  respiration  in  the  water  and  sediments  and  due  to 
diminished  saturation  dissolved  oxygen  concentration.  Conversely,  annual 
peaks  in  dissolved  oxygen  occur  in  winter  when  low  temperature  diminishes 
respiration  in  the  water  column,  light  stimulates  benthic  algal  production,  and 
saturation  concentration  is  maximum.  Our  model  indicates  minimum  daily- 
average  dissolved  oxygen  concentrations  of  2  to  3  gm  m  3  occur  in  upper 
Indian  River  (Figures  9-2,  9-3).  In  lower  Indian  River  and  Rehoboth  Bay,  we 
typically  indicate  minimum  daily- average  dissolved  oxygen  concentrations  of 
6  gm  m3  (Figures  9-5  to  9-9).  Observations  indicate  virtually  no  dissolved 
oxygen  concentrations  below  2  gm  m'3.  In  lower  Indian  River  and  Rehoboth 
Bay,  however,  a  substantial  proportion  of  observations  fall  below  our  6  gm  m  3 
minimum.  A  portion  of  the  discrepancy  lies  in  the  problem  of  comparing  grab 
samples  with  daily  average  model  output,  as  previously  described.  Several 
aspects  of  the  comparisons  remain  problematic,  however.  Observations  at 
most  stations  indicate  dissolved  oxygen  was  lower  in  1989  and  1990  than  in 
1988  (Figures  9-4  through  9-9).  We  cannot  explain  this  phenomenon;  our 
model  indicates  dissolved  oxygen  is  equivalent  in  the  three  summers.  Our 
model  also  agrees  well  with  the  summer  1988  observations.  Observations  at 
most  stations  indicate  minimum  dissolved  oxygen  in  1990  occurred  at  the 
beginning  of  May,  at  a  water  column  temperature  less  than  10  C°,  rather  than 
later  in  the  year.  We  believe  the  data  base  merits  careful  examination  before 
conclusions  are  drawn  regarding  the  ability  of  the  model  to  reproduce  the 
annual  cycle  of  dissolved  oxygen. 

A  mass-balance  feature  of  the  model  allows  computation  of  internal  mater¬ 
ial  sources  and  sinks.  These  are  compared  with  the  external  loads  quantified 
earlier  (Chapter  3).  On  an  annual  basis,  nonpoint  sources  are  the  largest  nitro¬ 
gen  loads  to  the  Indian  River/Rehoboth  Bay  system  (Figures  9-32  through 
9-34).  Atmospheric  loads  are  the  second  largest  source  while  point  sources  are 
the  least  of  the  sources  quantified.  Most  of  the  nitrogen  load  to  the  system  is 
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Figure  9-32.  Annual  Mass  Balance,  1988 


Chapter  9  Water  Quality  Modal  Appfication  and  Analysis 


1989 


NITROGEN 


-2610.41 


-200  -100 


-11864.41 


-73S2.9PZ 

l - - — l - 

-15000  -7500 


772765.3 

U94.1 


SEDIMENTS 
ATMOSPHERE 
POINT  SOURCE 
12667  8 NONPOINT  SOURCE 
OCEAN  EXCHANGE 


o 

kg/DAT 


PHOSPHORUS 


SEDIMENTS 


ATMOSPHERE 
POINT  SOURCE 
NONPOINT  SOURCE 
OCEAN  EXCHANGE 


0 

kg/DAY 


CARBON 


SEDIMENTS 
PRIMARY  PI 


19223  3  rruliARY  PRODUCTION 

J9223.3  pQINT  SQURfc 

7}  10047  4  nonpoint  source 

J  10047.4  QC  EXCHflN6E: 


0 

kg/OAY 


Figure  9-33.  Annual  Mass  Balance,  1989 
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Figure  9-34.  Annual  Mass  Balance,  1990 
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lost  through  the  mouth  to  the  ocean.  The  remainder  is  lost  to  the  sediments 
where  it  is  buried  or  denitrified  to  gaseous  form. 

During  years  of  dry  and  normal  hydrology,  the  largest  source  of  phosphorus 
to  Indian  River/Rehoboth  Bay  is  import  through  the  inlet  from  the  ocean  (Fig¬ 
ures  9-32,  9-34).  In  extremely  wet  years,  however,  nonpoint-  source  loads 
exceed  import  from  the  ocean  (Figure  9-33).  Point  source  loads  are  less  than 
the  oceanic  import  or  nonpoint  sources.  No  atmospheric  load  is  thought  to 
exist  The  sediments  represent  the  only  phosphorus  sink.  Since  no  process 
analogous  to  denitrification  exists  for  phosphorus,  the  ultimate  sink  is  burial  to 
deep,  inactive  sediments. 

During  years  of  dry  and  normal  hydrology,  the  largest  source  of  carbon  to 
Indian  River/Rehoboth  Bay  is  primary  production  by  algae  (Figures  9-32, 

9-34).  The  carbon  produced  by  algae  is  oxidized  in  the  water  column,  as 
BOD,  or  in  the  sediments,  as  sediment  oxygen  demand.  In  a  wet  year,  how¬ 
ever,  a  surprisingly  large  amount  of  carbon  enters  from  nonpoint  sources  so 
that  nonpoint-sources  comprise  the  largest  fraction  of  oxygen  demand  (Fig¬ 
ure  9-33).  Point-source  carbon  loads  are  negligibly  small.  Most  of  the  carbon 
load  to  the  system  settles  to  the  sediments  where  it  is  oxidized,  producing 
sediment  oxygen  demand,  or  buried  to  deep,  inactive  sediments.  A  lesser 
fraction  of  the  carbon  load  is  exported  through  the  mouth  to  the  ocean. 


Summer-Average  Water  Quality 

Spatial  representations  of  water  quality  (Figures  9-35  through  9-44)  were 
obtained  by  averaging  model  results  for  June,  July,  and  August  1990,  an  aver¬ 
age  hydrologic  year.  Summer  was  selected  for  spatial  representation  since  it  is 
the  season  in  which  water  quality  problems  most  often  occur. 

The  model  indicates  chlorophyll  (Figure  9-36)  is  maximum  near  the  Mills- 
boro  spillway.  High  chlorophyll  concentrations  also  occur  in  the  headwaters 
of  other  major  branches.  Chlorophyll  generally  declines  away  from  runoff 
sources  and  is  least  near  die  Indian  River  inlet  The  chlorophyll  distribution  is 
governed  by  die  distributions  of  ammonium  (Figure  9-38),  nitrate  (Fig¬ 
ure  9-39),  and  dissolved  phosphate  (Figure  9-42).  Nitrate  and  ammonium  are 
abundant  relative  to  the  half-saturation  concentration  for  nitrogen  uptake 
(Table  6-2),  throughout  most  of  Indian  River.  In  Rehoboth  Bay  and  Indian 
River  Bay,  however,  nitrate  and  ammonium  are  nearly  depleted.  The  spatial 
distribution  of  nitrate  reflects  the  source  in  runoff  and  wasteloading  and  the 
sink  in  the  ocean.  The  ammonium  distribution  indicates  a  major  source  is 
sediment  release.  Dissolved  phosphate  is  maximum  in  the  inlet,  indicating 
concentration  in  this  region  is  strongly  influenced  by  the  ocean  boundary  con¬ 
dition.  Phosphate  is  also  high  near  Millsboro  and  the  Rehoboth  Beach  outfall. 
Minimum  phosphate  concentrations  occur  midway  between  the  inlet  and  Mills¬ 
boro  and  in  portions  of  Rehoboth  Bay.  The  high  chlorophyll  levels  in  upper 
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Figure  9-36.  Modeled  Chlorophyll,  1990  Summer  Average 
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Figure  9-38.  Modeled  Ammonium,  1990  Summer  Average 
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Figure  9-41 .  Modeled  Total  Nitrogen,  1990  Summer  Average 
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Figure  9-43.  Modeled  Total  Phosphorus,  1990  Summer  Average 
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Indian  River  reflect  the  availability  of  nitrogen  and  phosphorus  for 
algal  growth.  Phosphorus  is  the  limiting  nutrient  in  this  region  (Figures  9-21 
through  9-23).  In  lower  Indian  Riyer  and  Rehoboth  Bay,  phosphorus  is  avail¬ 
able  to  support  algal  growth  but  nitrogen  is  limiting  (Figures  9-24  through 
9-28).  The  transition  from  phosphorus  limitation  in  headwaters  to  nitrogen 
limitation  near  the  ocean  interface  is  a  typical  estuarine  characteristic. 

The  distribution  of  total  nitrogen  (Figure  9-41)  shows  clearly  the  transition 
from  sources  along  the  shore  to  the  sink  in  the  ocean.  Total  phosphorus 
exhibits  die  same  trend  (Figure  9-43). 

Results  of  the  model  indicate  anoxia  is  not  a  permanent  feature  anywhere  in 
the  system  (Figure  9-44).  Summer-average  dissolved  oxygen  exceeds  6  gm  m 
3  everywhere.  Observations  and  the  model  indicate,  however,  that  concentra¬ 
tions  as  low  as  2  gm  m'3  occur.  Excursions  below  average  occur  on  a  diurnal 
basis  due  to  algal  respiration,  and  on  the  time  scale  of  days  due  to  flow  events 
and  similar  phenomenon.  Highest  dissolved  oxygen  concentrations  coincide 
with  highest  chlorophyll  concentrations  indicating  the  role  of  algal  production 
in  dissolved  oxygen  production.  Minimum  dissolved  oxygen  occurs  about 
halfway  between  MiUsboro  and  Indian  River  inlet.  The  origin  of  the  minimum 
is  not  apparent  but  chlorophyll  is  also  low  in  this  region  suggesting  the  mini¬ 
mum  is  due  to  limited  algal  production. 

A  substantial  and  significant  difference  exists  between  annual  mass  balances 
(Figures  9-32  through  9-34)  and  summer-average  mass  balances  (Figures  9-45 
through  9-47).  During  summers  of  dry  and  normal  hydrology  (Figures  9-45, 
9-47)  sediments  are  the  largest  nitrogen  source  to  the  system,  through  the 
mechanism  of  ammonium  release.  Even  in  an  extremely  wet  summer,  sedi¬ 
ments  rival  nonpoint  sources  of  nitrogen  (Figure  9-46).  The  summer  mass 
balance  contrasts  with  the  annual  balance  in  which  sediments  are  a  net  nitro¬ 
gen  sink.  A  similar  balance  exists  for  phosphorus.  The  sediment  release  of 
phosphorus  is  sufficient  to  convert  the  system  into  a  net  exporter  of  total  phos¬ 
phorus  during  the  summer.  By  contrast,  on  an  annual  basis,  the  sediments  are 
a  phosphorus  sink  and  the  system  imports  phosphorus  from  die  ocean.  The 
role  of  sediments  in  summer  in  Indian  River/Rehoboth  Bay  is  consistent  with 
the  role  of  sediments  in  larger,  deeper  estuaries.  The  net  nutrient  release  in 
summer  is  due  to  two  factors.  High  temperature  which  enhances  diagenetic 
(decay)  processes  in  the  sediments  and  seasonal  turbidity  increase  which 
shades  out  benthic  algae. 
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Conclusions 

Hydrodynamic  Model  Application 

The  CH3D  hydrodynamic  model  simulated  transport  in  Indian  River/ 
Rehoboth  Bay  for  three  years:  1988,  1989,  and  1990.  Two-dimensional  depth- 
average  transport  was  calculated  on  a  grid  containing  over  2000  cells  roughly 
50  m  by  100  m  in  lateral  and  longitudinal  extent  Integration  time  step  was 
30  seconds. 

Performance  of  die  model  was  verified  by  comparison  to  observed  tide 
records,  current  measurements,  and  salinity  data.  The  model  reproduced 
observed  surface  levels  throughout  the  system  with  an  accuracy  of  a  few  centi¬ 
meters.  Current  measures,  concentrated  in  the  passage  between  Rehoboth  Bay 
and  lower  Indian  River,  were  also  well-represented  by  the  model.  The  pre¬ 
dicted  spatial  salinity  distribution  compared  well  with  observations  collected  in 
three  independent  surveys. 

Computation  of  transport  in  Indian  River/Rehoboth  Bay  represents  the 
current  state  of  the  art  in  this  field.  The  transport  computations  are  entirely 
suited  for  use  in  a  eutrophication  model  of  the  system. 


Water  Quality  Model  Application 

The  water  quality  model  was  applied  in  a  continuous  mode  to  simulate 
conditions  in  Indian  River/Rehoboth  Bay  for  the  years  1988-1990.  Transport 
for  water  quality  simulation  was  provided  by  the  CH3D  hydrodynamic  model. 
Other  water-quality  forcing  functions  inducted  nonpoint-source  loading  of 
nutrients  and  organic  matter,  point-source  loading  of  nutrients  and  organic 
matter,  atmospheric  nutrient  loading,  sediment-water  nutrient  and  oxygen 
exchanges,  oceanic  material  exchanges,  and  seasonal  light  attenuation.  Quanti¬ 
fication  of  loads  and  other  forcing  functions  was  performed  specifically  for 
this  study  and  employed  die  most  recent  available  information.  The  model 
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provided  a  consistent  representation  of  water  quality  throughout  the  system 
throughout  the  year. 

The  water  quality  model  operated  on  a  tempo 'a  i  >cale  of  =300  seconds  and 
a  spatial  scale  of  *500  m.  Forcing  functions  were  not  resolved  on  these 
scales,  however.  Consequently,  comparisons  of  model  results  with  observa¬ 
tions  were  more  favorable  when  model  and  observations  were  temporally  and 
spatially  aggregated.  Performance  statistics  for  the  Indian  River/Rehoboth  Bay 
model  application  were  nearly  at  the  state  of  the  art  for  most  substances  simu¬ 
lated.  Unfavorable  statistics  noted  for  two  substances  were  attributed  to 
unique,  highly-localized  conditions  and  to  dubious  observations. 

The  model  package  is  entirely  suitable  for  employment  to  meet  the  objec¬ 
tives  stated  at  the  initiation  of  the  project.  Existing  point-source  and  nonpoint- 
source  loads  are  quantified  and  located  throughout  the  system.  These  can  be 
adjusted  to  develop  total  maximum  daily  loads  for  nutrients  and  organic  mat¬ 
ter.  Diurnal  dissolved  oxygen  variations  are  obtained  by  a  process  that  quanti¬ 
fies  diurnal  range  as  a  function  of  predicted  daily-average  dissolved  oxygen 
and  chlorophyll  concentrations. 


Recommendations 

Application  to  the  existing  model  package  to  the  existing  data  base  is  opti¬ 
mal.  No  improvement  of  model-data  agreement  is  possible  without  collection 
of  additional  data.  Representation  of  die  system  cannot  be  improved  without 
revisions  to  the  model  formulations.  Suggested  improvements  to  the  data  base 
and  the  model  follow. 


Data  Collection 

Water-column  observations  should  be  collected  at  consistent  locations  in 
Indian  River/Rehoboth  Bay  and  at  consistent  depths.  The  collection  of  surface 
samples  only  should  be  replaced  with  a  procedure  that  provides  samples  repre¬ 
sentative  of  tire  entire  water  column.  A  recommended  protocol  is  to  collect 
samples  at  mid-depth  when  depth  is  less  than  1  m.  Collect  surface  and  bottom 
samples  in  water  of  depth  between  1  and  2  m.  Collect  surface,  mid-depth,  and 
bottom  samples  in  water  of  depth  greater  than  2  m.  Salinity,  temperature,  and 
dissolved  oxygen  are  low-cost  measures,  often  observed  in  situ.  These  can  be 
measured  at  multiple  depths  easily.  If  resources  for  analysis  are  limited,  sam¬ 
ples  of  nutrients  and  organic  matter  can  be  composited  into  one  depth-average 
sample  prior  to  analysis. 

The  existing  data  base  is  highly  biased  towards  the  months  April  through 
September.  While  these  are  the  months  in  which  water  quality  problems  are 
most  severe,  lack  of  observations  in  the  remaining  months  limits  analyses  that 
might  relate  summer  water  quality  conditions  to  phenomenon  that  occur  in 
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colder  months,  for  example  load  events  or  winter  diatom  blooms.  Regular 
sampling  throughout  the  calendar  year  is  recommended. 

The  BODS  observations  currently  collected  are  no  longer  favored  as  a 
component  of  long-term  estuarine  data  bases.  Neither  is  BOD  a  state  variable 
in  state-of-the-art  eutrophication  models.  The  BOD  observation  should  be 
replaced  or  supplemented  with  observations  of  total  organic  carbon.  The  total 
organic  carbon  analyses  are  necessary  in  observations  of  point-source  and 
nonpoint-source  loads  as  well  as  the  water-column.  Direct  analysis  of  particu¬ 
late  mid  dissolved  organic  carbon  fractions  will  greatly  facilitate  future  model 
applications. 

The  detection  levels  of  ammonium,  nitrate,  and  phosphate  in  the  STORET 
data  base,  1988-1990,  are  unsuited  for  the  levels  of  these  nutrients  often  found 
in  Indian  River/Rehobo th  Bay.  A  large  fraction  of  the  observations  are  below 
detection  level.  Much  unproved  detection  levels  are  possible,  as  demonstrated 
in  the  College  of  Marine  Studies  data  base. 

Light  extinction  greatly  affects  water  quality  in  Indian  River/Rehoboth  Bay. 
The  existing  data  base,  however  is  assembled  from  observations  of  disk  visibil¬ 
ity  (secchi  depth).  The  conversion  of  disk  visibility  to  extinction  relies  on  an 
approximate  conversion.  Moreover,  we  found  the  disk  visibility  observations 
in  the  STORET  data  base  to  be  so  unreliable  we  could  not  use  them.  Direct, 
in-situ,  measurement  of  extinction  with  a  photometer  or  similar  instrument  is 
recommended. 

The  present  program  of  sediment-water  nutrient  and  oxygen  flux  measures 
should  be  continued.  Measures  should  be  collected  at  several  times  throughout 
the  year.  Concurrent  with  die  flux  measures,  water-column  light  extinction 
and  light  intensity  at  the  sediment-water  interface  should  be  recorded.  In  view 
of  the  role  of  benthic  algae,  chlorophyll  measurements  in  surficial  sediments 
should  be  added  to  the  sampling  program. 


Mechanistic  Sediment  Model 

One  component  of  the  Chesapeake  Bay  model  package  applied  to  Indian 
River/Rehoboth  Bay  was  a  mechanistic,  predictive  model  of  sediment-water 
nutrient  and  oxygen  exchanges.  Employment  of  the  benthic  sediment  model 
was  discontinued,  however,  since  it  did  not  account  for  activity  of  benthic 
algae  noted  in  Indian  River.  The  mechanistic  model  was  replaced  with  empiri¬ 
cal  models  calibrated  to  available  observations.  The  empirical  models 
performed  well  in  relating  existing  fluxes  to  existing  conditions,  primarily 
temperature  and  light.  Inherent  in  tire  empirical  models,  however,  was  the 
assumption  that  basic  fluxes  will  not  be  altered  as  a  result  of  management 
activities.  The  empirical  models  also  related  sediment-water  fluxes  to  ambient 
conditions  only.  The  models  did  not  relate  summer  nutrient  fluxes  to  loads  or 
benthic  algal  blooms  that  occurred  the  previous  winter.  Reactivation  of  the 
mechanistic  sediment  model  is  recommended.  In  order  to  employ  die  model, 
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however,  it  must  be  modified  to  account  for  the  effects  of  benthic  algae.  A 
reasonable  starting  point  for  the  modifications  is  adaptation  to  benthic  algae  of 
existing  kinetics  formulations  that  represent  water-column  phytoplankton. 


Multi-Layer  Hydrodynamic  Model 

The  data  base  for  Indian  River/Re  hobo  th  Bay  provided  little  information  on 
vertical  variations  in  water  quality.  Most  observations  were  surface  only. 
Examination  of  limited  available  observations,  collected  in  Indian  River  (Tyler 
1989),  indicated  vertical  variations  in  temperature,  salinity,  and  dissolved  oxy¬ 
gen  were  small,  typically  less  than  1  C°  for  temperature,  1.5  ppt  for  salinity, 
and  0.15  gm  m3  for  dissolved  oxygen.  No  distinct  pyncnocline  was  evident 
and  bottom-water  anoxia  was  absent.  Representation  of  water  quality  condi¬ 
tions  in  Indian  River/Rehobo  th  Bay  did  not  require  a  multi-layer  model. 

The  slight  salinity  stratification  suggests,  however,  that  classic  density- 
driven  circulation  exists,  to  some  extent,  in  die  reach  from  Millsboro  spillway 
to  the  Indian  River  inlet  The  density-driven  circulation  is  superimposed  on 
the  twice-daily  flow  reversals  induced  by  tides  and  is  characterized  by  net 
motion  downstream  near  the  surface  and  upstream  near  die  bottom.  The 
strength  of  die  circulation  depends  on  runoff  volume  and  salinity  at  the  inlet 
The  density-driven  circulation  is  subject  to  alteration  by  local  wind  conditions. 
At  present  effects  of  multi-layer  circulation,  to  the  extent  that  they  exist  are 
represented  by  a  dispersion  term  in  the  mass-transport  equatioa  This  repre¬ 
sentation  is  adequate,  as  evidenced  by  agreement  in  modeled  and  observed 
salinity  at  most  locations  and  times.  Still,  some  improvement  in  modeled 
circulation  may  be  expected  by  employment  of  a  multi-layer  model.  In  view 
of  the  shallow  depth  in  Indian  River/Re hoboth  Bay,  two  layers  would  suffice. 
No  multi-layer  modeling  can  be  conducted,  however,  until  an  extensive  data 
base  of  vertical  salinity  and  temperature  observations  is  assembled.  Long-term 
current  measures  at  several  locations  in  the  vertical  are  also  desirable. 
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model.  Other  water-quality  forcing  functions  inducted  nonpoint-source  loading  of  nutrients  and  organic  matter, 
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13.  Concluded 


Project  objectives  required  a  model  package  suitable  for  development  of  total  maximum  daily  loads,  point- 
source  waste  load  allocations,  and  nonpoint-source  load  allocations  of  nutrients  and  organic  substances.  Addi¬ 
tional  requirements  were  that  the  model  package  operate  in  a  continuous  multiyear  mode  and  provide  information 
on  diurnal  dissolved  oxygen  variations.  The  model  package  was  also  to  provide  an  organized  framework  for 
collection  and  employment  of  additional  observations  in  the  study  system.  The  model  package  described  in  this 
report  is  entirely  suitable  for  employment  to  meet  these  objectives. 


